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1.  Introduction 


This  report  describes  accomplishments  made  in  four  areas  of  research  related  to  robust, 
fixed-structure  control  system  design  and  analysis: 

1)  the  formulation  of  robust,  fixed-structure  control  design  problems  in  terms  of  a 
Riccati  equation  feasibility  and  the  development  of  probability-one  homotopy 
algorithms  for  its  solution  ; 

2)  a  comparison  of  descent  and  continuation  techniques  for  optimal  reduced-order 
control  design  and  an  investigation  of  the  best  bases  in  which  to  represent  the 
reduced-order  controller; 

3)  the  development  of  parallel  processing  techniques  to  implement  probability-one 
homotopy  algorithms  for  reduced-order  H2///00  control  design;  and 

4)  the  development  and  implementation  of  an  object-oriented  programming  approach 
for  the  implementation  of  interior  point  methods  to  solve  linear  matrix  inequalities 
(LMI’s). 

Below,  we  motivate  each  of  the  four  areas  of  research  and  briefly  describe  key 
accomplishments.  The  journal  paper  titles  are  used  as  headings  and  we  list  both  the 
journal  papers  that  have  been  submitted  and  the  corresponding  conference  papers.  Each 
of  the  journal  papers  is  included  in  an  Appendix. 


2.  Probability-One  Homotopy  Algorithms  for  Robust  Controller 
Synthesis  with  Fixed-Structure  Multipliers 

Journal  Paper 

E.  G.  Collins,  Jr.,  W.  M.  Haddad,  L.  T.  Watson,  and  D.  Sadhukhan,  “Probability-One 
Homotopy  Algorithms  for  Robust  Controller  Synthesis  with  Fixed-Structure  Multipliers,” 
International  Journal  of  Robust  and  Nonlinear  Control,  accepted  for  publication. 

Conference  Paper 

E.  G.  Collins,  Jr.,  W.  M.  Haddad,  and  L.  T.  Watson,  “Fixed-Architecture,  Robust  Control 
Design  Using  Fixed-Structure  Multipliers,”  Proceedings  of  the  International  Federation 
of  Automatic  Control,  to  appear  in  June  1996. 

During  the  past  two  decades,  major  advancements  have  been  made  in  robust  control 
theory.  Building  upon  theory,  the  structured  singular  value  (SSV)  was  defined  as  a 
nonconservative  robustness  measure  for  the  analysis  of  linear  systems  with  arbitrary  phase, 
multiple-block  uncertainty.  The  supremum  of  the  structured  singular  value  over 


nonnegative  frequencies  is  the  inverse  of  the  multivariable  stability  margin.  The  initial 
developments  in  structured  singular  value  theory  focussed  on  uncertainty  with  arbritrary 
phase  (often  called  “complex  uncertainty”)  and  hence,  although  less  conservative  than 
theory,  could  still  yield  very  conservative  robustness  bounds  for  systems  with  parametric 
uncertainty.  This  led  to  the  development  of  mixed  (i.e.,  real  and  complex)  structured 
singular  value  (MSSV)  theory  which  considers  block-diagonal  uncertainty  with  both 
complex  and  real  scalar  parametric  elements. 

Parallel  research  addressed  the  issue  of  real  parameter  uncertainty  using  absolute 
stability  theory  such  as  Popov  analysis  and  was  developed  by  recognizing  the  relationship 
between  sector  bounded  nonlinearities  and  interval  bounds  on  linear  uncertainties.  This 
work  was  soon  seen  to  provide  an  upper  bound  for  the  MSSV.  In  fact,  in  contrast  to  the 
initial  work  on  the  MSSV,  this  research  provided  the  first  fixed-structure  multiplier 
versions  of  MSSV  theory.  A  unique  contribution  of  some  of  this  work  is  that  it  led  to  the 
development  of  upper  bounds  on  an  H2  cost  functional  over  the  uncertainty  set  under 
consideration.  By  optimizing  this  upper  bound  and  using  a  Riccati  equation  constraint, 
continuation  algorithms  have  been  developed  for  MSSV  controller  synthesis.  Note  that 
the  H2  approach  allows  the  direct  design  of  fixed-architecture  (e.g.,  reduced-order  or 
decentralized)  controllers  and  the  simultaneous  optimization  of  the  controller  and  (fixed- 
structure)  multipliers,  hence  avoiding  M-K  (i.e,  multiplier-controller)  iteration  schemes. 
However,  to  date  the  synthesis  algorithms  have  been  formulated  only  for  the  case  of  the 
the  Popov  multiplier.  In  addition,  the  algorithms  rely  on  an  initialization  scheme,  have  not 
used  the  prediction  capabilities  obtained  by  computing  the  Jacobian  of  the  homotopy  (or 
continuation)  map,  and  have  assumed  that  the  homotopy  curve  is  monotonic. 

A  similar  line  of  research  has  been  developed  independently  by  Safonov,  et.  al.  This 
work  also  provides  a  fixed-structure  multiplier  version  of  the  MSSV  but,  unlike  the 
approach  described  above,  this  approach  develops  multipliers  for  strictly  linear 
uncertainties.  The  associated  robustness  analysis  was  originally  formulated  in  terms  of  a 
convex,  frequency-domain  optimization  problem  but  has  recently  been  reformulated  in 
terms  of  a  (convex)  linear-matrix-inequality  (LMI)  problem.  These  results  have  led  to  the 
recognition  that  robust  control  design  can  be  approached  via  solving  a  (nonconvex) 
“bilinear  matrix  inequality”  (BMI).  This  approach  allows  the  design  of  fixed-architecture 
controllers  and  can  be  implemented  without  using  M-K  iteration.  To  obtain  a  reasonably 
sized  BMI,  the  multiplier  set  must  be  restricted  to  lie  in  the  span  of  a  stable  basis. 
However,  the  choice  of  this  basis  is  unclear  and  can  potentially  introduce  a  high  degree  of 
conservatism.  If  the  less  conservative  LMI  formulation,  requiring  the  use  of  unstable 
multipliers,  is  used,  the  resultant  BMI  is  of  very  high  dimension  due  to  the  introduction  of 
a  Lyapunov  inequality  of  the  dimension  of  the  closed-loop  system  to  ensure  closed-loop 
stability.  In  contrast,  the  robustness  analysis  results  using  a  Riccati  equation  formulation 
easily  extend  to  robust  control  design  without  placing  any  basis  restrictions  on  the 
multipliers  or  introducing  high  dimensionality. 

This  research  used  a  Riccati  equation  constraint  to  formulate  fixed-architecture,  robust 
control  design  methods  that  use  general  forms  of  the  fixed-structure  multipliers.  The 
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proposed  method  relies  on  the  development  of  an  artificial  cost  function.  This  cost 
function  also  includes  barrier  functions  to  enforce  positive  definite  constraints  (e.g.,  on  the 
Riccati  solution  P)  which  allows  the  constrained  optimization  problem  (the  constraints 
including  P>0)  to  be  converted  into  an  unconstrained  optimization  problem.  The  cost 
function  is  not  physically  meaningful  so  we  do  not  encounter  the  normal  problems 
associated  with  making  the  barrier  functions  small  as  the  last  step  of  the  optimization 
process.  If  the  barrier  terms  are  ignored  and  a  certain  term  is  added  to  the  cost  function, 
the  cost  function  becomes  an  H2  upper  bound. 

Due  to  the  positive  definite  constraint  on  the  Riccati  solution,  it  is  not  possible  to 
approach  the  solution  to  the  optimization  problem  using  standard  descent  methods. 
Hence,  we  develop  probability-one  homotopy  algorithms  to  find  the  solution.  This  class 
of  homotopy  algorithms  is  distinct  from  classical  continuation  algorithms  in  that  they 
follow  the  zero  curve  using  the  arc  length  parameter  and  not  the  actual  homotopy 
parameter  X.  This  allows  the  zero  curve  to  be  nonmonotonic  in  X  and  provides  additional 
numerical  robustness.  In  addition,  the  algorithms  developed  can  be  initialized  with  any 
stabilizing  compensator  and  admissible  multiplier,  in  contrast  to  previous  algorithms. 


3.  A  Comparison  of  Descent  and  Continuation  Algorithms  for  Ht 
Optimal,  Reduced-Order  Control  Design 

Journal  Paper 

E.  G.  Collins,  Jr.  and  D.  Sadhukhan,  “A  Comparison  of  Descent  and  Continuation 
Algorithms  for  H2  Optimal  Reduced-Order  Control  Design,”  submitted  to  the 
International  Journal  of  Control. 

Conference  Paper 

E.  G.  Collins,  Jr.  and  D.  Sadhukhan,  “A  Comparison  of  Descent  and  Continuation 
Algorithms  for  H2  Optimal  Reduced-Order  Control  Design,”  to  be  submitted  to  the  1997 
American  Control  Conference. 

One  of  the  deficiencies  of  modern  control  laws,  developed  by  simply  solving  a  pair  of 
decoupled  Riccati  equations,  in  particular,  linear-quadratic-guassian  (LQG)  control  and 
full-order  control,  is  that  the  resultant  control  laws  are  always  of  the  order  of  the 
design  plant.  These  techniques,  though  relatively  easy  to  implement  computationally,  do 
not  allow  the  designer  to  constrain  the  architecture  (e.g.  order  and  degree  of 
centralization)  of  the  controller.  Such  constraints  are  often  necessary  in  engineering 
practice  due  to  throughput  limitations  of  the  control  processors.  Reduced-order  control  is 
therefore  of  paramount  importance  in  practical  control  design.  This  research  focuses  on 
the  design  of  H2  optimal,  reduced-order  controllers. 

Two  main  approaches  have  been  developed  to  solve  the  H2  optimal,  reduced-order 
design  problem.  The  first  approach  attempts  to  develop  approximations  to  the  optimal 
reduced-order  controller  by  reducing  the  dimension  of  an  LQG  controller.  These  methods 


are  attractive  because  they  require  relatively  little  computation  and  should  be  used  if 
possible.  Unfortunately,  they  tend  to  yield  controllers  that  either  destabilize  the  system  or 
have  poor  performance  as  the  requested  controller  dimension  is  decreased  or  the 
requested  control  authority  level  is  increased.  Hence,  if  used  in  isolation,  these  methods 
do  not  yield  a  reliable  methodology  for  reduced-order  design.  In  addition,  these  methods 
do  not  extend  to  the  design  of  decentralized  controllers.  However,  it  should  be  mentioned 
that,  in  regards  to  reduced-order  control  design,  the  indirect  approaches  at  worst  are 
valuable  in  providing  good  initial  conditions  for  the  direct  approaches  described  below. 

In  contrast  to  controller  reduction,  direct  approaches  attempt  to  directly  synthesize  an 
optimal,  reduced-order  (or  decentralized)  controller  by  a  numerical  optimization  scheme. 
There  are  two  main  classes  of  parameter  optimization  approaches  to  direct  control  design. 
The  first  class  relies  on  the  use  of  descent  methods.  Algorithms  in  this  class  reduce  the 
Hi  cost  at  each  iteration.  The  second  class  relies  on  the  use  of  continuation  methods.  In 
contrast  to  the  descent  methods,  the  Hi  cost  is  not  necessarily  reduced  at  each  iteration.  It 
should  be  mentioned  that  continuation  algorithms  have  also  been  developed  to  solve  the 
“optimal  projection  equations,”  a  set  of  four  coupled  Lyapunov  and  Riccati  equations  that 
characterize  the  Hi  optimal,  reduced-order  compensator.  However,  this  approach  was  not 
be  considered  in  this  research. 

From  a  practical  design  perspective  it  is  important  to  determine  which  class  of  methods 
tends  to  be  more  numerically  robust.  As  with  the  vast  majority  of  numerical  methods  for 
nonconvex  optimization  problems,  answers  to  these  questions  are  extremely  difficult  to 
prove  analytically.  Instead,  we  must  rely  on  numerical  experimentation  to  observe  trends. 
Hence,  in  this  research  the  behavior  of  some  standard  descent  methods  (i.e.,  steepest 
descent,  conjugate  gradient,  BFGS  Quasi-Newton,  and  Newton's  method)  were  compared 
to  the  corresponding  behavior  of  the  continuation  algorithm  of  by  considering  design  for 
three  reduced-order  control  design  problems  appearing  in  the  literature.  The  results 
clearly  indicate  that  the  continuation  algorithm  tends  to  be  more  numerically  robust  and  is 
most  efficient  when  the  controller  is  constrained  to  a  tridiagonal  form. 


4.  Cost-Effective  Parallel  Processing  for  H2/H00  Controller  Synthesis 

Journal  Paper 

Y.  Ge.,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “Cost-Effective  Parallel  Processing  for  HilH^ 
Controller  Synthesis,”  submitted  to  the  International  Journal  of  Control. 

Conference  Paper 

Y.  Ge.,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  "Distributed  Homotopy  Algorithms  for 
H2/H„  Controller  Synthesis,"  in  Parallel  Processing  for  Scienticfic  Computing,  P.  E. 
Bjorstand,  J.  R.  Gilbert,  M.  V.  Mascagni,  R.  S.  Schreiber,  H.  D.  Simon,  V.  J.  Tqrczor, 
and  L.  T.  Watson,  eds.,  SIAM,  Philadephia,  PA,  pp.  84-89,  1995. 


H2/H„  mixed-norm  controller  synthesis  is  an  important  and  interesting  technique  in 
modem  control  design  which  provides  the  means  for  simultaneously  addressing  H2  and 
performance  objectives.  In  practice  such  controllers  provide  both  nominal  performance 
(via  suboptimal  H2)  and  robust  stability  (via  H^).  Hence  mixed-norm  synthesis  provides  a 
technique  for  trading  off  performance  and  robustness,  a  fundamental  objective  in  control 
design. 

The  H2/H0.  mixed-norm  problem  has  been  addressed  in  a  variety  of  settings.  One 
treatment  utilized  an  H2  cost  bound  as  the  basis  for  an  auxiliaiy  nonconvex  constrained 
minimization  problem,  which  is  veiy  difficult  to  solve  without  the  global  convergence  of 
homotopy  methods.  A  successful  homotopy  algorithm  based  on  the  Ly  form 
parametrization  has  previously  been  developed. 

The  H2/H^  control  design  algorithms  can  be  used  for  controller  design  of  systems  such 
as  the  benchmark  four  disk  system.  This  system  is  especially  representative  of  the  type  of 
vibration  control  problems  that  arise  in  industrial  problems  involving  rotating 
turbomachinery.  H2/H„  design  can  be  used  to  develop  controllers  that  are  robust  with 
respect  to  unmodeled  dynamics  and  also  guarantee  a  certain  measure  of  nominal 
performance. 

It  should  be  mentioned  that  H2/H„  theory  has  been  used  previously  to  develop  complex 
structured  singular  value  (CSSV)  synthesis  formalisms  that  a  priori  fix  the  structure  of 
both  the  D-scales  and  the  controller.  Hence,  an  extension  of  the  algorithms  here  will 
enable  fixed-structure  CSSV  controller  synthesis  that  blends  H2  and  /L,  performance 
objectives. 

Practical  applications  often  lead  to  large  dense  systems  of  nonlinear  equations  which 
are  time-consuming  to  solve  on  a  serial  computer.  For  these  systems,  parallel  processing 
may  be  the  only  feasible  means  to  achieving  solution  algorithms  with  acceptable  speed. 
One  economical  way  of  achieving  parallelism  is  to  utilize  the  aggregate  power  of  a 
network  of  heterogeneous  serial  computers.  In  industrial  environments  where  interactive 
design  is  often  the  practice,  the  parallel  code  can  be  easily  incorporated  into  interactive 
software  such  as  MATLAB  or  Mathematica  with  proper  setup  of  the  network  computers. 
To  the  engineering  users  the  design  environment  is  identical.  However,  the  computations 
are  faster. 

The  most  expensive  part  of  the  previously  developed  H2IH0.  homotopy  algorithm  is  the 
computation  of  the  Jacobian  matrix,  which  can  be  parallelized  easily  to  run  across  an 
Ethernet  network  with  little  modification  of  the  original  sequential  code,  and.  which  has 
relatively  large  task  granularity.  There  is  a  trade-off  between  the  programming  effort  and 
the  speedup  of  the  parallel  program.  To  obtain  a  better  speedup,  other  parts  of  the 
homotopy  algorithm,  such  as  finding  the  solution  to  the  Riccati  equations  and  the  QR 
factorization  to  compute  the  kernel  of  the  Jacobian  matrix,  need  to  be  parallelized  as  well. 
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In  this  sutdy  the  homotopy  algorithm  for  controller  synthesis  was  parallelized  to 

run  on  a  network  of  workstations  using  PVM  (Parallel  Virtual  Machine)  and  on  an  Intel 
Paragon  parallel  computer,  under  the  philosophy  that  as  few  changes  as  possible  are  to  be 
made  to  the  sequential  code  while  achieving  an  acceptable  speedup.  The  parallelized 
computation  is  that  of  the  Jacobian  matrix,  which  is  carried  out  in  the  master-slave 
paradigm  by  functional  parallelism;  that  is,  each  machine  computes  a  different  column  of 
the  Jacobian  matrix  with  its  own  data.  Unless  the  Riccati  equation  solver  is  parallelized, 
there  is  a  large  amount  of  data  needed  for  each  slave  process  at  each  step  of  the  homotopy 
algorithm.  To  avoid  sending  too  many  large  messages  across  the  network  or  among 
different  nodes  on  the  Intel  Paragon,  all  slave  processes  repeat  part  of  the  computation 
done  by  the  master  process,  which  therefore  decreases  the  speedup  of  the  parallel 
computation. 

The  speedups  of  the  parallel  code  were  compared  as  the  number  of  workstations 
increases  or  as  the  number  of  nodes  increases  on  an  Intel  Paragon  and  as  the  size  of  the 
problem  varies.  A  reasonable  speedup  can  be  achieved  using  an  existing  network  of 
workstations  compared  to  that  of  using  an  expensive  parallel  machine,  the  Intel  Paragon. 
It  was  demonstrated  that  for  a  large  problem,  the  approach  of  using  a  network  of 
workstations  to  parallelism  is  feasible  and  practical,  and  provides  an  efficient  and 
economical  computational  method  to  parallelize  a  homotopy  based  algorithm  for  H2/H^ 
controller  synthesis  in  a  workstation-based  interactive  design  environment. 


5.  An  Object-oriented  Approach  to  Semidefinite  Programming 

Journal  Paper 

Y.  Ge,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “An  Object-oriented  Approach  to 
Semidefinite  Programming,”  submitted  to  the  International  Journal  of  Computer 
Mathematics. 

Conference  Paper 

Y.  Ge,  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “An  Object-oriented  Approach  to 
Semidefinite  Programming,”  to  be  submitted  to  the  1997 American  Control  Conference. 

Object-oriented  design  and  programming  has  been  a  major  theme  in  software 
engineering  in  recent  years.  Traditional  design  or  classical  design,  which  has  been  the 
main  software  design  paradigm  until  about  the  mid  80’ s,  concentrates  on  the  actions  that  a 
system  has  to  take  and  decomposes  the  sytem  into  separte  units  or  modules  according  to 
their  functionalities.  In  object-oriented  design  a  system  to  be  modeled  is  viewed  as  a 
collection  of  objects,  each  of  which  has  its  own  attributes  and  the  operations  performed  on 
an  object  or  functions  acting  on  an  object  are  also  defined  in  one  syntactic  unit.  Objects 
communicate  by  passing  messages  or  by  calling  functions  from  other  objects  which 
provide  services.  Object-oriented  design  is  developing  an  object-oriented  model  of  a 
system  and  can  be  realized  (implemented)  by  object-oriented  programming  using 
languages  such  as  C-I-+,  FORTRAN  90,  or  Smalltalk. 


The  advantages  of  object-oriented  design  and  programming  have  been  described 
widely  elsewhere.  A  short  summary  is  provided  here.  First,  an  object  is  an  independent 
entity  that  is  encapsulated  in  one  syntactic  unit.  The  definition  of  an  object  consists  of  the 
definition  of  the  attributes  of  the  object  along  with  operations  that  can  be  performed  on 
the  object  and  the  services  or  function  calls  provided  by  the  object.  Encapsulation  hides 
the  implementation  details  of  an  object  and  makes  the  program  easier  to  modify.  Any 
subsequent  changes  to  the  program  can  be  localized,  making  the  resulting  program  more 
easily  maintained. 

The  second  advantage  is  information  hiding.  Definitions  of  an  object  which  need  not 
be  known  to  other  objects  are  inaccessible  to  other  objects,  preventing  them  from  being 
changed  accidentally.  In  other  works,  information  hiding  makes  implementation  details  of 
an  object  inaccessible  to  other  objects.  However,  the  designer  has  the  freedom  to  decide 
what  to  hide  and  what  not  to  hide. 

The  third  advantage  is  code  reuse.  Inheritance  enables  the  definition  of  a  new  object, 
which  can  be  viewed  as  a  subclass  of  an  existing  object,  without  having  to  repeat  some  of 
the  details.  The  new  object  can  inherit  attributes  or  operations  from  its  ancestor. 
Inheritance  is  one  way  to  support  reuse  of  existing  objects.  There  are  different  kinds  of 
reuse  in  object-oriented  programming;  inheritance  is  only  one  of  them. 

One  of  the  most  popular  object-oriented  programming  languages  is  C-I-+,  which  was 
used  to  implement  the  algorithm  of  this  research.  Some  of  the  reasons  why  C+-t-  is  so 
widely  used  are  upward  compatibility  with  C,  design  emphasis  on  efficiency  and 
performance,  and  the  availability  of  many  useful  libraries  and  tools.  For  instance,  the  Gnu 
C-l-t-  compiler  and  other  tools  are  available  on  a  wide  range  of  platforms  and  provide  good 
performance,  programming  environments,  and  reasonable  compliance  with  ANSI 
standards. 

There  are  many  available  libraries  such  as  IML++,  SparseLib-H-,  STL,  and  others 
which  emphasize  numerical  computation.  One  notable  package  is  LAPACK-h-,  developed 
by  Dongarra  et  al.,  which  is  a  C-t--!-  interface  to  LAPACK  and  BLAS.  It  has  been  shown 
that  performance  of  programs  using  the  package  is  comparable  to  calling  LAPACK  and 
BLAS  directly,  anc  can  at  the  same  time  reap  benefits  of  object-oriented  programming. 

This  research  developed  and  analyzed  object-oriented  design  and  implementation  of  an 
algorithm  for  semidefmite  programming.  Semidefinite  programming  refers  to  minimizing  a 
linear  function  subject  to  a  linear  matrix  inequality.  That  is 

minimize  x 


are  symmetric  matrices.  F(x)>0  means  that  F(x)is 


ceR’^\  andFQ,...,F^ 
positive  semidefmite. 

Many  problems  in  controls  engineering  can  be  cast  in  terms  of  a  semidefmite 
programming  problem.  Since  a  semidefmite  programming  problem  is  a  convex 
optimization  problem,  which  can  be  solved  by  interior  point  methods,  it  has  attracted  the 
attention  of  many  researchers  in  interior  point  methods.  There  is  a  C  implementation  of  a 
primal-dual  algorithm  for  the  semidefinite  problem.  A  C++  implementation  of  that  primal- 
dual  algorithm  was  developed  in  this  resarch  to  explore  the  possible  benefits  of  object- 
oriented  design.  Because  of  the  similarity  of  the  primal-dual  algorithm  with  other  interior 
point  algorithms  for  solving  the  semidefmite  programming  problem,  the  design  and 
implementation  methodology  developed  in  this  research  can  be  easily  modified  and  applied 
to  other  interior  point  algorithms. 

The  performance  of  a  C++  implementation  of  the  primal-dual  algorithms  for 
semidefmite  programming  is  compared  with  the  existing  C  implementation.  While  the 
CPU  times  of  the  two  implementations  are  comparable  to  each  other,  the  C++  version 
offers  the  advantages  mentioned  earlier  in  this  section.  Segments  of  the  code  are  used  to 
illustrate  object-oriented  features  of  the  implementation. 
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Summary 


Continuation  algorithms  that  avoid  multiplier-controller  iteration  have  been  developed  earlier 
for  fixed-architecture,  mixed  structured  singular  value  controller  synthesis.  These  algorithms  have 
only  been  formulated  for  the  special  case  of  Popov  multipliers  and  rely  on  an  ad  hoc  initialization 
scheme.  In  addition,  the  algorithms  have  not  used  the  prediction  capabilities  obtained  by  computing 
the  Jacobian  matrix  of  the  continuation  (or  homotopy)  map,  and  have  assumed  that  the  homotopy 
zero  curve  is  monotonic.  This  paper  develops  probability-one  homotopy  algorithms  based  on  the 
use  of  general  fixed-structure  multipliers.  These  algorithms  can  be  initialized  using  an  arbitrary 
(admissible)  multiplier  and  a  stabilizing  compensator.  In  addition,  as  with  all  probability-one  al¬ 
gorithms,  the  homotopy  zero  curve  is  not  assumed  to  be  monotonic  and  prediction  is  accomplished 
by  using  the  homotopy  Jacobian  matrix.  This  approach  also  appears  to  have  some  advantages  over 
the  bilinear  matrix  inequality  (BMI)  approaches  resulting  from  extensions  of  the  LMI  framework 
for  robustness  analysis. 

Keywords:  Fixed-structure  multipliers,  multivariable  stability  margin,  LMFs,  BMPs,  homotopy 
algorithms 

Running  Title:  Robust  Controller  Synthesis  with  Fixed-Structure  Multipliers 


1.  Introduction 


During  the  past  two  decades,  major  advancements  have  been  made  in  robust  control  theory. 
Building  upon  theory,  the  structured  singular  value  (SSV)^’  ^  was  defined  as  a  nonconservative 
robustness  measure  for  the  analysis  of  linear  systems  with  dynamic,  arbitrary  phase,  multiple- 
block  uncertainty.  The  supremum  of  the  structured  singular  value  over  nonnegative  frequencies  is 
the  inverse  of  the  multivariable  stability  margin  (see  References  .3,  4  and  the  references  therein). 
The  initial  developments  in  structured  singular  value  theory  focussed  on  dynamic  uncertainty  with 
arbritrary  phase  (often  called  “complex  uncertainty”)  and  hence,  although  less  conservative  than 
theory,  could  still  yield  very  conservative  robustness  bounds  for  systems  with  parametric  un¬ 
certainty.  This  led  to  the  development  of  mixed  (i.e.,  real  and  complex)  structured  singular  value 
(MSSV)  theory^’  ®  which  considers  block-diagonal  uncertainty  with  both  dynamic  and  real  scalar 
parametric  elements. 

Parallel  research  addressed  the  issue  of  real  parameter  uncertainty  using  absolute  stability  the¬ 
ory  such  as  Popov  analysis^ and  was  developed  by  recognizing  the  relationship  betw’een  sector 
bounded  nonlinearities  and  interval  bounds  on  linear  uncertainties.  This  work  was  soon  seen  to 
provide  an  upper  bound  for  the  MSSV.®’  In  fact,  in  contrast  to  the  initial  work  on  the  MSSV, 
this  research  provided  the  first  fixed-structure  multiplier  versions  of  MSSV  theory.  A  unique  con¬ 
tribution  of  some  of  this  work  is  that  it  led  to  the  development  of  upper  bounds  on  an  H2  cost 
functional  over  the  uncertainty  set  under  consideration.  By  optimizing  this  upper  bound  and 
using  a  Riccati  equation  constraint,  continuation  algorithms  have  been  developed  for  MSSV  con¬ 
troller  synthesis.^'* A  related  algorithm  for  complex  stuctured  singular  value  (CSSV)  controller 
synthesis  is  given  in  Reference  17.  Note  that  the  H2  approach  allows  the  direct  design  of  fixed- 
architecture  (e.g.,  reduced-order  or  decentralized)  controllers  and  the  simultaneous  optimization 
of  the  controller  and  (fixed-structure)  multipliers,  hence  avoiding  M-K  (i.e.,  multiplier-controller) 
iteration  schemes.  However,  to  date  the  synthesis  algorithms  have  been  formulated  only  for  the  case 
of  the  the  Popov  multiplier.  In  addition,  the  algorithms  rely  on  an  ad  hoc  initialization  scheme,  have 
not  used  the  prediction  capabilities  obtained  by  computing  the  Jacobian  matrix  of  the  homotopy 
(or  continuation)  map,  and  have  assumed  that  the  homotopy  curve  is  monotonic. 

A  similar  line  of  research  has  been  developed  independently  in  Reference  18-20.  This  work  also 
provides  a  fixed-structure  multiplier  version  of  the  MSSV  but,  unlike  tlie  approach  described  in  Ref¬ 
erences  7-10,12,  this  approach  develops  multipliers  for  strictly  linear  uncertainties.  The  associated 
robustness  analysis  was  originally  formulated  in  terms  of  a  convex,  frequency-domam  optimiza¬ 
tion  problem  but  has  recently  been  reformulated  in  terms  of  a  (convex)  linear-matrix-inequality 


(LMI)  problem. These  results  have  led  to  the  recognition  that  robust  control  design  can  be  ap¬ 
proached  via  solving  a(nonconvex)  “bilinear  matrix  inequality”  This  approach  allows 

the  design  of  fixed-architecture  controllers  and  can  be  implemented  without  using  M-K  iteration. 
To  obtain  a  reasonably  sized  BML  the  multiplier  set  must  be  restricted  to  lie  in  the  span  of  a  stable 
basis. ““  However,  the  choice  of  this  basis  is  unclear  and  can  potentially  introduce  a  high  degree 
of  conservatism.  If  the  less  conservative  LMI  formulation,  requiring  the  use  of  unstable  multipli¬ 
ers,  is  used,  the  resultant  BMI  is  of  very  high  dimension  due  to  the  introduction  of  a  Lyapunov 
inequality  of  the  dimension  of  the  closed-loop  system  to  ensure  closed-loop  stability. In  contrast, 
the  robustness  analysis  results  using  a  Riccati  equation  formulation  easily  extend  to  robust  control 
design  without  placing  any  basis  restrictions  on  the  multipliers  or  introducing  high  dimensionality. 

This  paper  uses  a  Riccati  equation  constraint  to  formulate  fixed-architecture,  robust  control 
design  methods  that  use  general  forms  of  the  fixed-structure  multipliers.  The  proposed  method 
relies  on  the  development  of  an  artificial  cost  function.  This  cost  function  also  includes  barrier 
functions  to  enforce  positive  definite  constraints  (e.g.,  on  the  Riccati  solution  P)  which  allows 
the  constrained  optimization  problem  (the  constraints  including  P  >  0)  to  be  converted  into  an 
unconstrained  optimization  problem.  The  cost  function  is  not  physically  meaningful  so  we  do  not 
encounter  the  normal  problems  associated  with  making  the  barrier  functions  small  at  the  last  step 
of  the  optimization  process.  (See  Reference  25  for  a  discussion  of  this  negative  feature  of  standard 
barrier  function  methods.  )  If  the  barrier  terms  are  ignored  and  a  certain  term  is  added  to  the  cost 
function,  the  cost  function  becomes  an  H2  upper  bound. 

Due  to  the  positive  definite  constraint  on  the  Riccati  solution,  it  is  not  possible  to  approach 
the  solution  to  the  optimization  problem  using  standard  descent  methods.  Hence,  we  develop 
probability-one  homotopy  algorithms'^’  to  find  the  solution.  This  class  of  homotopy  algorithms 
is  distinct  from  classical  continuation  algorithms'^  in  that  they  follow  the  zero  curve  using  the 
arc  length  parameter  and  not  the  actual  homotopy  parameter  A,  This  allows  the  zero  curve  to 
be  nonmonotonic  in  A  and  provides  additional  numerical  robustness.  In  addition,  the  algorithms 
developed  here  can  be  initialized  with  any  stabilizing  compensator  and  admissible  multiplier,  in 
contrast  to  the  algorithms  of  References  14-17. 

1,1.  Notation  and  Definitions 

Let  TZ  and  C  denote  the  real  and  complex  numbers,  and  the  real  and  complex 

m  X  m  matrices,  let  (•)^  and  (•)*  denote  transpose  and  complex  conjugate  transpose,  let  “Re”  and 
“Im”  denote  real  and  imaginary  part,  and  let  or  I  denote  the  nxn  identity  matrix.  Furthermore 


we  write  crjnax( *)  for  the  maximum  singular  value,  “tr"  for  the  trace  operator,  and  A/  >  0  (M  >  0)  to 
denote  the  fact  that  the  Hermitian  matrix  M  is  nonnegative  (positive)  definite.  The  Hermitian  and 
skew-Hermitian  parts  of  an  arbitrary  complex  square  matrix  G  are  defined  by  He  G  =  \(G  +  <?*) 
and  Sh  G  =  respectively.  Finally,  vec(-)  denotes  the  standard  column  stacking  operator. 

Next,  we  establish  certain  key  definitions  used  later  in  the  paper.  Let  n(s)  and  d{s)  be  poly¬ 
nomials  in  s  with  real  coefficients.  A  function  g(s)  of  the  form  g(s)  —  n(s)/d{s)  is  called  a  real 
rational  function.  The  function  ^^(.s)  is  called  proper  (resp.,  strictly  proper)  if  deg  n(.s)  <  deg  d{s) 
(resp.,  deg  n{s)  <  deg  d{s)),  where  “deg"  denotes  the  degree  of  the  respective  polynomials.  A  real- 
rational  matrix  function  is  a  matrix  whose  elements  are  rational  functions  with  real  coefficients. 
Furthermore,  a  transfer  function  G(  s)  is  called  proper  (  resp.,  strictly  proper)  if  every  element  of  Gis) 
is  proper  (resp.,  strictly  proper).  In  this  paper  we  assume  all  transfer  functions  are  real-rational 
matrix  functions.  Also,  we  define  G*(s)  =  ^^(-*.5)  for  transfer  functions  G(s), 

All  asymptotically  stable  transfer  function  is  a  transfer  function  each  of  whose  poles  is  in  the  open 
left  half  plane.  Finally,  a  Lyapunov  stable  transfer  function  is  a  transfer  function  each  of  whose  poles 
is  in  the  closed  left  half  plane  with  semi-simple  poles  on  the  imaginary  axis.  Let  G{s}  ~ 
denote  a  state  space  realization  of  a  transfer  function  G(s),  that  is,  G(s)  =  C{sl  —  A)~^B  +  D, 

A  square  transfer  function  G(s)  is  called  positive  real  (resp.,  generalized  positive  real)^^  if  G(s) 
is  Lyapunov  stable  and  He  G{s)  is  nonnegative  definite  for  Re [.s]  >  0  (resp.,  G(s)  has  no  imaginary 
poles  and  He  G(ju;)  is  nonnegative  definite  for  all  o:  G  7^).  A  square  transfer  function  is  called 
strictly  positive  real^^  (resp.,  strictly  generalized  positive  real)  if  G{s)  is  asymptotically  stable  and 
He  G{ju))  is  positive  definite  for  u;  G  7^  (resp.,  G{s}  has  no  imaginary  poles  and  He  Gifu)  is 
positive  definite  for  a;  G  Tv).  A  square  transfer  function  G{s)  is  strongly  positive  real  (resp.,  strongly 
generalized  positive  real)  if  it  is  strictly  positive  real  (resp.,  strictly  generalized  positive  real)  and 
>  0,  where  D  =  G{oo).  Note  that  although  a  minimal  realization  of  a  positive  real  transfer 
function  is  stable  in  the  sense  of  Lyapunov,  a  minimal  realization  of  a  generalized  positive  real 
transfer  function  may  be  unstable. 

1.2.  Paper  Organization 

Section  2  presents  the  general  framework  for  robustness  analysis  with  fixed-structure  multipliers 
and  briefly  discusses  the  BMI  approach  to  fixed- architecture,  robust  control  synthesis.  Section  3 
demonstrates  that  an  alternative  formulation  to  robust  synthesis  is  in  terms  of  a  Riccati  equation 
feasibility  problem.  Section  4  develops  probability-one  homotopy  algorithms  to  solve  the  Riccati 
equation  feasibility  problem.  Section  5  specializes  the  results  to  the  special  case  of  the  Popov 
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multiplier  and  presents  a  numerical  example.  Finally,  Section  6  presents  conclusions  and  directions 
for  further  research. 


2.  Multiplier  Methods  in  Robust  Analysis 


In  this  section  we  review  the  framework  for  mixed  uncertainty  robustness  analysis  with  fixed- 
structure  multipliers.  The  exposition  generally  follows  that  presented  in  References  12,19-21.  VVe 
begin  by  considering  the  standard  uncertainty  feedback  configuration  of  Figure  1.  where  G'(s)  € 


^mxm  jg  asymptotically  stable  and  G(s)  ^ 
belongs  to  the  set 


It  is  assumed  that  the  uncertainty  A 


V 

=  {A  =  block-diag(  Ai, ...,  Ap) ;  A,  €  T,-,  <Tmax(  A,)  <  7,  i  =  1,  ....p,  ki  =  m},  (2.1) 

i=l 

where  J,  C  denotes  the  internal  structure  of  the  uncertainty  block  A,  and  7  >  0.  For 

example  J,-  may  be  given  by  any  of  the  following  five  sets: 


(2.2) 

A”  A  {A,-  =  6ilk,  :  6i  €  C}, 

(2.3) 

A«i  A  |A.-  =  8ih,  :  8i  €  7^}, 

(2.4) 

=  {Ai  €71^'^'^'  :  A,  =  Af}, 

(2.5) 

{A,  =  [  ~1‘' 

(2.6) 

Note  that  A^,  A^^,  and  A^^^  are  standard  in  the  literature,  corresponding  respectively  to  com¬ 
plex  matrix  block  uncertainty,  repeated  complex  scalar  uncertainty,  and  repeated  real  scalar  un¬ 
certainty.  A*'^  is  symmetric,  real  matrix  block  uncertainty,  while  can  be  used  to  describe 
uncertainty  in  the  imaginary  and  real  parts  of  a  structural  system  represented  in  real  normal  form. 
If  the  uncertainty  is  of  the  form  described  by  A^'^  or  A^,  it  is  possible  to  represent  the  uncertainty 
by  A^^^.  However,  as  discussed  in  Reference  12  this  reformulation  leads  to  increased  conservatism 
and  numerical  complexity.  The  ensuing  discussion  is  not  restricted  to  these  forms  of  uncertainty,  but 
they  are  important  special  cases  and  wiU  be  used  to  provide  concrete  illustrations  of  the  subsequent 
concepts. 


To  state  the  multivariable  absolute  stability  criterion  for  A  €  we  define  the  sets  of  Hermitian. 
frequency-dependent,  scaling  matrix  functions  by 

Vi  =  {Di  :  i7^  U  oo  —  :  Diijuj)  >  0, 

Di(ju>)Ai  =  AiDiiju).  LJ  €  72,  A,  €  I;,  i  =  1,  (2.7) 

Mi  =  {Ni  :  j72  U  00  -  =  Ni*{j^), 

Ni{ju})Ai  =  ArNiiJu).  u;  G  72,  A,-  €  lu  i  =  1, -.p)-  (2.8) 

Furthermore,  define  the  sets  Mi  and  M  of  multiplier  transfer  functions  by 

Mi  =  {Miis)  =  Di{s)  -h  Qi(.s) :  Diiju)  €  Vi,  Qi(jui)  =  juNiijuj),  Niijuj)  €  Mi,  i  =  1, 

(2.9) 

M  =  {M(s)  e  :  M{s)  =  block-diag(M,(.s). ...,Mp{-s)),  Mi(s)  e  Mi,  i  =  (2.10) 

Note  that  in  (2.9),  Di(joj)  =  He  Miiju)  >  0  and  Ni(ju)  =  -jSli  Mi{juj).  Furthermore.  M(s)  e  M 


satisfies 


and  is  not  necessarily  stable. 


He  M{ju)  >0,  w  €  72  U  oo. 


(2.11) 


Theorem  1.  Suppose  that  G'(.s)[/  -  7G(.s)]  ^  is  asymptotically  stable.  If  there  exists  M{s)  6 


M  such  that 


where 


Re[M{juj)T^{j(jj)]  >  0.  u;  €  72U  oo. 


T^is)^[I  +  ^G{s)][I-'rGis)]-\ 


(2.12) 


(2.13) 


then  the  negative  feedback  interconnection  of  G(s)  and  A  is  asymptotically  stable  (or,  equivalently, 
det(/  -I-  G(j(jj)A)  0,  w  €  72)  for  all  A  €  A-^. 

Proof.  A  rigorous  proof  of  this  result  is  given  in  Reference  12.  Similar  results  are  presented  in 
References  19-21.  A 


Remark  1.  Note  that^^ 


A  +  7S(/-7jD)-^C 


(2.14) 
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Using  the  coprime  factorization  result  presented  in  Reference  31,  it  follows  that  M(s)  can  be 
factored  as 

A/{.s)  =  [A/B(^r'MA(.s),  (2.15) 

where  both  Ma(-s)  and  Mb(s)  are  asymptotically  stable  and  nonunique.  In  practice,  stable  Ma(s) 
and  -1/b(  s)  satisfying  (2.15)  can  be  computed  using  the  approach  pioneered  in  Reference  .32  which  is 
also  given  in  Reference  31.  In  Reference  20  (2.15)  is  used  to  prove  the  following  important  corollary 
to  Theorem  1. 

Corollary  1.  Assume  iV/(.s)  G  M,  and  A/a(-5)  aad  Mb{s)  are  asymptotically  stable  transfer 
function  matrices  satisfying  (2.15).  Then  (2.12)  in  Theorem  1  holds  if  and  only  if 


Ee[MA(ju:)T^(j(^)MB(joj)]  >  0,  w  e  IZU  oo. 


(2.16) 


Corollary  1  allows  us  to  develop  robust  controllers  based  on  positive  real  theory.  Such  a  method 
is  developed  in  Section  3. 

We  now  characterize  Mi  for  T,-  equal  to  the  sets  defined  by  (2.2)-(2.6).  These  characterizations 
are  restatements  of  results  given  in  Reference  12.  Multiplier  sets  corresponding  to  A^^,  and 
are  given  in  References  19,  20,  The  set  corresponding  to  differs  from  that  given  here. 

The  following  characterizations  are  useful  in  constructing  state  space  realizations  of  the  multipliers. 
In  particular,  for  I,-  =  A^ 

Mi  =  :  Re[mj(ju;)]  >  0,  lm[mi(ju)]  =  0,  a;  G  7^  U  oo}  ,  (2.17) 

while  for  =  A^^ 

Mi  =  {m.(5)  €  :  He  M.-(jw)  >  0,  M{jw)  =  M*(jw),  w  €  U  00}  .  (2.18) 

Note  that  if  we  denote  Mi(s)  =  k-)-,  M{juj)  =  M*(  ju}  implies 

Im[mjj(iw)]  =  0,  mjkijw)  =  mkji-jw),  j  7^  k. 

Furthermore,  for  I,:  =  A^^* 

Mi  =  {a/,(s)  €  :  He  M(jw)  >  0,  w  G  72  U  00}  ,  (2.19) 

while  for  J,-  =  A^ 


Mi  =  {m,(s)Ijt;  ;  mi{s)  G  C,  Re[mi(ja;)]  >  0,  u;  G  72  U  00} 


(2.20) 


7 


Finally,  for  J,  = 

Mi  =  {Mi(s)  =  DU)  +  Q{.?)  :  D{s)  = 


dn(s)  di2{s) 
-duis)  dii(s) 


Qi^^) 


.  He  D(ju})  >  0, 


9ll(-s)  9i2(s) 

912(5)  -911(5) 

Re  Q(ju>)  =  0,lm[f(ii(jw)]  =  Re[rfi2(iti;)]  =  0,  u;  6  72.U  00}. 


(2.21) 


2.1.  The  Structured  Singular  Value  and  Robust  Performance 

For  a  multiple  block-structured  uncertainty  set  I.  with  possibly  repeated  scalar  elements,  com¬ 
plex  scalar  elements,  real  blocks,  and  complex  blocks. the  structured  singular  value  of  a  complex 
matrix  G{ju!)  is  defined  by 


l.i{G{juj))  =  inf 

\Ael 

where  by  convention  /i(G'(ya;))  =  0  if  there  does  not  exist  A  £  X  such  that  det(/  -f  G{jui)A)  =  0. 
The  structured  singular  value  nonconservatively  characterizes  the  robust  stability  of  the  uncertainty 
feedback  system  of  Figure  1,  as  stated  by  the  following  theorem.^’ 


{o-maxCA)  :  det(/  +  G{ju)A)  =  0} 


Theorem  2.  Suppose  G(s)  is  asymptotically  stable.  Then  the  negative  feedback  interconnec¬ 
tion  of  G(s)  and  A  is  asymptotically  stable  for  all  A  €  Ay  if  and  only  if 

^{G(jijj))  <  7~^,  u;  G  7^  U  00.  (2.22) 

Remark  2.  The  parameter  7m  =  1/sup  fi{G(ju))  is  the  multivariable  stability  margin.^ 

i-/>0 

Next  define 

/tabs(G(jt<;) )  =  inf{7  >  0  :  there  exists  M{-)  €  M  such  that  He[ili‘(yw)Ty(ya;)]  >  0,  w  €  7^  U  00}. 

(2.23) 

Then  the  following  holds. 

Theorem  3.  For  w  G  7^  U  00,  let  G{juj)  be  a  complex  matrix.  Then 

H{G{ju))  <  n^hs{G{joj)). 

Proof  A  rigorous  proof  is  given  in  Reference  12.  Similar  results  are  stated  in  References  19,  20. 


A 


The  significance  of  the  above  theorem  is  that  it  allows  us  to  consider  both  robust  stability  and 
robust  performance  in  the  same  setting.*^  Hence,  as  was  proved  for  structured  singular  value  analysis 
with  purely  complex  uncertainty,’  robust  performance  can  be  ensured  by  appropriate  inclusion  of 
a  “fictitious"  full  complex  uncertainty  block. 

2.2.  Linear  Matrix  Inequality  and  Riccati  Equation  Characterizations  of  (Gen¬ 
eralized)  Positive  Real  Matrix  Functions 


Theorem  1  characterizes  robustness  in  terms  of  the  strictly  generalized  positive  real  condition 
(2.12)  while  Corollary  1  relies  on  the  strictly  positive  real  condition  (2.16).  Hence,  to  implement  the 
robustness  test  of  Theorem  1  using  state  space  computations  requires  state  space  characterizations 
of  strictly  generalized  positive  real  and  strictly  positive  real  transfer  functions. 


State  space  conditions  for  strictly  positive  real  transfer  functions  are  given  in  Reference  .33,  but 
®  include  observability  and  rank  conditions  which  are  difficult  to  incorporate  into  numerical  schemes. 

Hence,  the  lemma  below  provides  state  space  characterizations  of  strongly  generalized  positive 
real  matrices  and  strongly  positive  real  matrices,  special  cases  respectively  of  strictly  generalized 
positive  real  matices  and  strictly  positive  real  matrices. 


Lemma  1.  Let  G'(-?)  be  square  with  G(s) 
foUowing  statements  are  equivalent: 


'  A 

B  ' 

C 

D 

where  (A,B)  is  controllable.  Then  the 


1.  G'(.s)  is  strongly  generalized  positive  real. 


2.  There  exist  c  >  0  and  symmetric  P  such  that 


-A^P-PA  -PB  +  C'^ 

-B^ P  +  C  (D-  €l)  +  {D-  elf 


>  0. 


Furthermore,  if  (1)  and  (2)  hold,  then: 


4. 

•5. 


(A,C)  is  observable  if  and  only  if  P  is  nonsingular. 

G(5)  is  strongly  positive  real  if  and  only  if  A  is  asymptotically  stable  and  P  >  0.  In  this  ca,se, 
{A.C)  is  observable  if  and  only  if  P  >  0. 


6. 


If  (A,  C)  is  observable,  G(s)  is  strongly  positive  real  if  and  only  if  Z>  +  >  0  and  there  exist 

P  >  0  and  c  >  0  such  that 


0  =  A^P  +.  PA  +  (P^P  -  cf{D  +  D^)  \P^P  -  C)  +  el. 


(2.24) 
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Proof.  The  equivalence  of  statements  1  and  2  is  shown  in  Reference  19.  Statements  3  and  4 
are  proved  in  Reference  29  and  statement  5  is  proved  in  Reference  9.  A 

Lemma  1  shows  that  strongly  generalized  positive  reality  and  generalized  positive  reality  is 
equivalent  to  the  existence  of  a  certain  solution  to  an  LML  Furthermore,  the  latter  condition  (under 
the  assumption  that  (yl,  B.  C\  D)  is  a  minimal  state  space  realization)  is  equivalent  to  the  existence 
of  a  certain  solution  to  a  Riccati  equation.  This  Riccati  equation  is  subsequently  used  to  develop 
probability-one  homotopy  maps  for  robust  controller  synthesis  with  fixed-structure  multipliers. 

2.3.  Bilinear  Matrix  Inequality  Approaches  to  Robust  Controller  Synthesis 

Both  References  19,  21  describe  ways  of  using  the  LMTs  corresponding  to  generalized  positive 
real  tests  to  develop  robustness  tests  expressed  in  terms  of  the  existence  of  a  solution  to  an  LML 
The  key  is  to  represent  the  multiplier  A/(.s)  in  the  strictly  generalized  positive  real  test  (2.12)  of 
Theorem  1  in  such  a  way  that  its  parameters  appear  linearly  in  the  corresponding  state  space  test 
of  Lemma  1.  When  the  LMI  robustness  analysis  results  are  generalized  to  fixed-architecture,  robust 
control  synthesis,  a  BMI  results. 

A  straightforward  way  of  achieving  the  desired  representation  of  M{s)  is  described  in  Reference 
21  and  is  based  on  earlier  ideas  given  in  Reference  20.  This  approach  requires  expressing  the 
multiplier  M(s)  in  terms  of  a  basis  expansion.  In  particular, 

r 

M{.s)  =  Y,i5jM^^\,s),  (2.25) 

j=i 

where  i3j  >  0  and  M^s)  6  M,  j  =  1, r.  A  substantial  weakness  of  this  approach  is  the  diflRicultv 
of  choosing  the  basis  ....  )}.  Specifically,  for  a  given  multiplier  M(s)  €  M  and  a 

given  basis,  the  approach  in  Reference  21  does  not  guarantee  that  there  exists  nonnegative  /3j  such 
that  M($)  is  given  by  (2.25). 

The  approach  proposed  in  Reference  19  does  not  suffer  from  these  weaknesses.  It  is  based  on 
the  following  result. 

Lemma  2.  Equation  (2.12)  is  satisfied  for  some  transfer  matrix  €  Af  if  and  only 

if  there  exist  a  real  polynomial  matrix  Af(.s)  6  M  for  which  (2.12)  holds.  Furthermore,  if  Af  (“^(s) 
is  factored  cis  Af(°)(s)  =  j(o)(^A(®i(s)  where  is  a  real  polynomial  matrix  and  d(°^(s)  is 

a  scalar- valued  real  polynomial,  then  the  degree  of  M(s)  need  not  be  greater  than  the  sum  of 
degrees  of  and  d(®i(a),  and  in  fact  one  can  choose  M(s)  =  d^°)(— 5)A^®)(s).  In  addition,  if 

M(s)  is  of  order  2n,  then  for  all  n*^^-order  real  polynomials  d(s)  having  no  zeros  on  the  joj  axis, 


M{s)  =  €  M  and  satisfies  (2,12). 


Remark  3.  Lemma  2  allows  one  to  restrict  the  multiplier  search  to  27^^^  order,  real  polynomial 
matrices  To  obtain  state  space  realizable  transfer  functions  we  can  consider 

M{s) 


d(s)d(~.sy 

where  d{.s)  is  an  arbitrary  n‘^-order  polynomial  having  no  zeros  on  the  jw-axis. 


(2.26) 


Remark  4.  Note  that  since  the  zeros  of  d(-s)  are  the  mirror  images  of  the  zeros  of  d{s)  about 
the  imaginary  axis,  M(s)  given  by  (2.26)  is  always  an  unstable  multiplier. 

This  latter  approach  is  powerful  but  as  eluded  to  in  Remark  4  always  produces  an  unstable 
multiplier.  When  extended  to  fixed  architecture,  robust  control  design,  this  approach  results  in 
a  BMI  with  very  high  dimension  since  it  requires  the  introduction  of  a  Lyapunov  inequality  of 
dimension  equal  to  that  of  the  closed-loop  system  to  ensure  closed-loop  stability.  We  hence  begin 
the  development  of  an  alternative  scheme  without  these  limitations. 

3.  A  Riccati  Equation  Approach  to  Robust  Controller  Synthesis  with  Fixed- 
Structure  Multipliers 

We  now  give  exclusive  attention  to  robustness  tests  that  are  expressed  in  terms  of  positive  real 
conditions,  as  opposed  to  generalized  positive  real  conditions.  Hence,  we  will  consider  robustness 
tests  corresponding  to  Corollary  1.  The  developments  here  differ  dramatically  from  those  in  the 
previous  section  since  we  use  the  Riccati  equation  characterization  of  strong  positive  reality  given 
in  Lemma  1  instead  of  the  LMI  characterization  of  strong  generalized  positive  reality  which  is  also 
given  in  Lemma  1. 

\\‘e  focus  on  the  uncertainty  structures  corresponding  to  the  sets  and  A*^^,  defined  re¬ 
spectively  by  (2.2)  and  (2.4).  Hence,  we  develop  fixed-structure  multiplier  tests  for  the  complex, 
block-structured  uncertainty  considered  by  classical  complex  structured  singular  value  analysis  ^ 
and  the  real,  diagonal  uncertainty  considered  by  classical  real  structured  singular  value  analysis.®’  ® 
Although  not  detailed  here,  the  uncertainty  structures  corresponding  to  the  sets  A^^,  A^'^,  and  A'^, 
may  also  be  considered  in  the  framework  of  this  section.  In  addition,  the  results  may  be  extended 
in  a  straightforward  manner  to  mixed  uncertainty  sets. 

The  resulting  numerical  algorithms  are  probability-one  homotopy  algorithms.  A  detailed  dis¬ 
cussion  of  the  numerical  algorithms  is  reserved  for  the  next  section.  Here,  we  lay  the  necessary 


foundation  for  these  algorithms.  We  begin  by  developing  a  constructive  characterization  of  multi¬ 
plier  factors  A/a(-s)  and  Msis)  corresponding  to  complex,  block-structured  uncertainty  and  real, 
diagonal  uncertainty.  These  constructions  are  essential  to  the  subsequent  developments. 


3.1.  Complex,  Block-Structured  Uncertainty 


From  (2.10)  and  (2.1/  )  it  follows  that  a  multiplier  corresponding  to  complex  block-diagonal 
uncertainty  is  given  by 


• 

Mis)  =  block-diag  {m.iis)Ik^,...,mpis)Ik^  , 

(3.1) 

Re  miiju)  >  0,  w  6  7^  U  00, 

(3.2) 

Im  niiijuj)  =  0.  a?  G  U  00. 

(3.3) 

•  Note  that  we  have  replaced  the  weak  inequality  in  (2.17)  with  a  strict  inequality  in  (3.2). 


Lemma  3.  For  uncertainty  structures  corresponding  to  (2.2)  there  exists  M{s)  with  no  zeros 
or  poles  on  the  imaginary  axis  and  satisfying  the  compatibility  conditions  (3.1)-(3.3)  and  the 
robustness  test  (2.12),  if  and  only  if  there  exists  Mis)  satisfying  (2.12)  and  (3.1)  such  that 

Mis)  =  M{-s)M(s),  (3.4) 


where  Mis)  has  no  poles  or  zeros  in  the  closed  right  half  plane. 


Proof.  Assume  there  exist  Mis)  satisfying  the  conditions  of  Lemma  3.  Then,  from  Lemma  2 
it  follows  that  there  exists  a  real  polynomial  Mis)  satisfying  (3.1)-(3.3)  and  (2.12).  Furthermore, 
the  restriction  (3.3)  also  requires  that  m,(5)  have  only  even  powers  of  s,  such  that  for  some  integer 
n,  niiis)  =  ■  This  implies  that 

n  Ti 

m,(s)  =  min  =  '^in 

k=l  k—l 

where  bk  =  here  denotes  the  square  root  in  the  right  half  plane. 

For  s  -V  oo  (3.2)  implies  that  for  some  real  positive  scalar  c,  (— l)”m,„  =  c^,  which  in  turn 
implies  that 

n 

m,(s)  =  =  ^t(-«)«t(5),  (3.6) 

fc=i 

where  n,(s)  has  no  zeros  in  the  closed  right  half  plane.  Note  that  if  Im(6jfe)  ^  0  there  exists  j  ^  k 
such  that  bj  =  bk.  It  follows  that  ra,(s)  =  +  bk)  is  a  real  polynomial.  Now  it  follows  from 
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(3.6)  that  M(s)  ==  N{—s)N{s)  where  N{s)  is  a  real  polynomial  matrix  with  no  zeros  in  the  closed 
right  half  plane.  If  d(s)  is  any  order,  real  polynomial  with  no  zeros  in  the  closed  right  half  plane 
and  M'(s)  ^  M(—s)M{s)  wdiere  M{s)  =  then  i\7'(.s)  has  no  zeros  or  poles  on  the  imaginary 
axis,  satisfies  (3.1)-(3.3)  and  from  Lemma  2  also  satisfies  (2.12). 

Finally,  note  that  each  M(s)  given  by  (3.1)  and  (3.4),  where  il7(.s)  has  no  poles  or  zeros  in  the 
closed  right  half  plane,  satisfies  (3.2)  and  (3.3).  A 

Remark  5.  Equation  (3.4)  corresponds  to  the  stable  coprime  factorization  of  A/(<s)  given  by 
(2.15)  with  A/a(-s)  =  M{s)  and  Mb{s)  =  M~^(s).  Furthermore,  if  M(s)  is  strictly  proper  then 
M{s)  given  by  (3.4)  is  strongly  positive  real  and  M{s)  and  both  have  state  space  realiza¬ 

tions. 


Remark  6,  M(s)  is  precisely  an  asymptotically  stable,  minimum  phase  transfer  function  rep¬ 
resentation  of  the  classical  D-scales  from  complex  structured  singular  value  analysis.^’  ^ 


From  Remark  5,  it  follows  that  if  we  denote  the  state  space  realization  of  strictly  proper  A/a(*s  ) 
in  (2.15)  by 


Ma(s} 


■  A 

B  ' 

c 

D 

(3.7) 


then,  we  can  simply  choose  A'/b(*5)  =  A/a  ^(5)  and  hence  using  a  standard  state  space  realization 
inversion  formula 


Msis) 


'  A-BD~^C 

BD~^  ' 

-D~'C 

D~' 

(3.8) 


Details  on  how  to  choose  (A.B.C.D)  to  enforce  the  block-diagonal  structure  of  A/a(-?)  are  given 
in  Reference  36. 

Now,  if  we  let  T^(s)  have  the  state  space  realization  (A^,  5^,  D^)  as  in  (2,14),  then,  referring 

to  (2.16), 

MA(s)T^is)MB(:S) 


A^ 

where 


— 


A-BD~^C  0  0 

0 

-BD^D~'^C  BB^  A 


,  B-y  = 


=  [  -DDyD~^C  DD^  C  ]  , 


BD~^ 

ByD~^ 

BD^D~' 


by  =  DD^D~^. 


(3.9) 

(3.10) 


The  next  theorem  which  considers  complex,  block-structured  uncertainty,  follows  immediately  from 
Corollary  1  and  Lemma  1. 


Theorem  4.  Let  Xi  =  i  =  1, _ and  suppose  G(.s)  is  asymptotically  stable.  If  there 

exist  (.4.  jB,  C,  jD),  P  >  0,  and  e  >  0  such  that  >  0  and 

0  =  +  PA^  +  P  -  C^)^{D^  +  D-,^)~^(b/p  -(:%)  +  e/; 

where  By,C\,  D-y  are  given  by  (3.9)  and  (3.10),  then  the  negative  feedback  interconnection  of 
6'(.s)  and  A  is  asymptotically  stable  for  all  A  €  A-,.. 

3.2.  Real,  Diagonal  Uncertainty 

Recall  from  (2.10)  and  (2.19)  that  a  multiplier  corresponding  to  real,  diagonal  uncertainty  (with 
possible  repeated  elements)  is  given  by 

Mis)  -  block-diag(iV/i(.s), ...,  A/p(s)), 

where 

He[A'/(7w)]  >  0.  ti;  €  PU  oo.  (3.11) 

Now.  if  we  let  Mis)  =  Afe*!?)”*  AfA(s)  as  in  (2.15),  then  it  follows  from  Corollary  1  that  if  we 
consider  only  strict  inequality  in  (3.11),  then  we  can  replace  (3.11)  with 

He[AfA(jw)A/B(ia/’)]  >  0.  w  €  72U  oo.  (3.12) 

Next  let  the  state  space  realizations  of  A'/a(s)  and  A/b(s)  be  denoted,  respectively,  by 


Ma(s) 


■  4a 

Bk  ■ 

Ca 

Dk  . 

,  A/b(s) 


Bb  ' 

Cb 

Db 

If  w'e  let  Tyis)  have  the  state  space  realization  (  j4.y,  B^,  D^)  as  in  (2.14),  then  referring  to  (2.16) 

MKis)Tyis)MBis) 

where 


^7,1 

B-i,\ 

^7,1 

0 

0 

CQ 

Bb 

47,1  — 

B^Cb  4.y  0 

1  ^7.1  ~ 

B^Db 

B\D^Cb  BaC^  Aa 

BaD^Db 

^7.1  = 

[  DaD^Db  DaC~^  Ca 

1  »  -^7,1  • 

=  DaB^Db- 

Similarly,  referring  to  (3.12), 


A/a(«)A/b(5) 


47,2 

By, 2 

^7,2 

By, 2 
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where 


Now.  let 


Ab  0 

•  5-/, 2  — 

Bb 

BaC'b  a  a 

BaDb 

[  BaCb  Ca 

•  -^'7,2  * 

=  DaDb- 

=  block-diagjyi.^.,!,  ^7  =  Wock-diag{jB.y,i, 

C-y  =  block-diag{C-y4,C-^,2},  =  block-diag{£>.y4,  .D-^,2}- 


#  The  next  theorem,  follows  immediately  from  Corollary  1  and  Lemma  1. 


(3.13) 

(3.14) 


Theorem  5.  Let  I,-  =  i  =  l,...,p,  and  suppose  G{s)  is  asymptotically  stable.  If  there 


p:.!" 


exist  (/Ia,5a,Ca,  L>a),  (Ab,  Bb,Cb,  Db),  P  >  0-  and  e  >  0  such  that  Dy  +  £>.,  >  0  and 

0  =  a/p  +  PA^  +  (B^'^p  -  (b^  +  bJf\B/p  -  C^)  +  €l. 


where  A^,  B^,C-f,  b-y  are  given  by  (3.13)  and  (3.14),  then  the  negative  feedback  interconnection  of 
G{s)  and  A  is  a^symptotically  stable  for  all  A  6  Ay. 

3.3.  Problem  Formulation 


Both  Theorems  4  and  5  provide  robustness  tests  in  terms  of  the  following  feasibility  problem. 

Riccati  Equation  Feasibility  Problem  (REFP).  Find  6  €  7^’,  c  >  0,  and  P  €  such 
that 


0  =  Ay'^{0)P  +  PAyie)  +  {By'^{e)P  -  Cy{9)f(by(e)  +  by^{e))  \B^{0)P-Cy{e))+€l,  (3.15) 


P>0,  by{0)  +  b'^{0)  >  0, 


(3.16) 


where  the  dimension  q  is  determined  by  the  multiplier  and  r  is  determined  by  both  the  multiplier 
and  the  nominal  plant  size.  In  Theorems  4  and  5,  0  corresponds  to  the  free  parameters  of  the 
matrices  providing  a  state-space  representation  of  the  multiplier  factors  M\(s)  and  Mb{s).  For 
0  example,  considering  Theorem  4,  if  all  of  the  elements  of  the  matrices  A,  B,  C,  and  D  in  (3.7)  and 

fj% 

(3.8)  are  free,  then  6  is  defined  by  6  =  ^vec^(.4),  vec^(R).vec^(C),  vec^(.D)j  . 

If  we  are  considering  control  design  for  a  plant  {Ap,Bp,Cp,Dp)  under  a  feedback  controller 
{Ac.  BcjCc),  then,  assuming  negative  feedback,  A  in  (2.14)  is  given  by 


A  = 


Ap  -BpCc 

BcCp  Ac  -  BcDpCc 


(3.17) 
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Hence,  in  Theorems  4  and  5  is  linear  in  the  controller  matrices.  The  controller  matrices  essen¬ 
tially  provide  extra  degrees  of  freedom  to  satisfy  the  Riccati  equation  constraint  (3.15  ).  To  illustrate, 
if  all  of  the  elements  of  the  matrices  A,  B.  C,  and  D  are  free,  and  all  of  the  matrices  Ac,  Cc 
are  also  free  then  9  is  defined  by  9  =  (vec^( A), vec^(5), vec^(C), vec^(.D), vec^(Ac), vec^(jBc)i 
vec^(/?c))^-  Note  that  Ay(9)^  By(9)^  C\{9),  and  D^{9)  are  generally  nonlinear  functions  of  9. 
Hence  it  is  not  possible  to  convert  the  REFP  to  an  L^II  feasibility  problem. 

We  approach  the  development  of  a  solution  technique  by  defining  the  following  artificial  cost 
function 

J(0,c,P)  =  a9'^9  +  ae^  +  tr  [D^(9)  +  i)^{9)]  +  Vp  tr  K  (3.18) 

where  a,  Vp,  and  are  positive  scalars.  This  barrier  cost  function  has  the  nice  properties  that 

1.  .7(0, e.P)  becomes  unbounded  if  one  of  the  eigenvalues  of  P  or  D^(9)  +  D^(9)  approaches 
zero  or  e  approaches  zero;  and 

2.  the  second  derivative  of  the  first  two  terms  with  respect  to  the  vector  [9,€]^  is  2aL 

Property  (1)  is  used  to  enforce  the  constraints  P  >  0,  D^.(9)  +  D^{9)  >  0,  and  c  >  0.  Property  (2) 
is  used  to  improve  the  numerical  conditioning  of  the  Hessian  matrix,  which  sometimes  improves 
the  numerical  robustness  of  the  associated  numerical  algorithms. 

Remark  7.  Replace  el  in  (3.15)  with  eR  where  R  =  E  and  E  is  the  output  matrix  corre¬ 
sponding  to  the  performance  variable  and  let  V  =  where  is  the  input  matrix  correspond¬ 

ing  to  the  plant  disturbance.  Then,  if  we  choose  J(9.e,P)  =  tr  PV  +  an  additional  nonnegative, 
multiplier-dependent  term,  then  ./(•)  is  an  upper  bound  on  an  H2  cost  functional.  This  type  of  cost 
functional  is  the  basis  for  the  continuation  and  homotopy  algorithms  of  References  14-17,  34,  37. 

Let  S  denote  the  set  of  solutions  (^,  c,  P)  to  the  REFP  and  consider  the  optimization  problem 

min  J(^,c,P)  subject  to  (3.15).  (3.19) 

We  will  attempt  to  find  {0,e,P)  €  by  solving  the  optimization  problem  (3.19).  That  is,  we  will 
search  for  (O.e,  P)  £S  that  is  an  extremal  to  (3.19). 

Note  that  we  are  assuming  that  if  S  is  nonempty,  (3.19)  has  a  solution.  However,  this  has  not 
been  proven.  A  sufficient  condition  to  guarantee  a  solution  is  that  the  feasible  level  set 

n(/3)  =  {(0,e,P)65:J(0,6,P)<^}, 
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is  compact  for  0  >  0.  Clearly,  n(/?)  is  bounded  due  to  the  terms  aO^O  +  ae^  and  the  positivity  of 
each  of  the  terms  in  (3.18).  However,  whether  n(/3)  is  closed  depends  on  the  closure  of  the  set  S. 
The  characterization  of  S  is  a  subject  of  future  research. 

To  characterize  the  extremals,  we  define  the  Lagrangian 


C(0,  €,P,Q)=^  J{0,  €.  P)  +  tr  QW{9,  c,  F), 


(3.20) 


where  WiO.e.P)  denotes  the  right  hand  side  of  (3.15)  and  Q  is  the  symmetric  Lagrange  multiplier. 
Note  that  the  constraints  (3.16)  are  absorbed  into  J  as  barriers.  The  necessary  conditions  for  a 
solution  to  (3.19)  are  given  by^® 

Hr  Hr 

(3.21) 


.  dC  ^  dc 
dO'  °  ■  c)c  ’ 


S/1  -  T  ~  ~T  -  T  ^ 

0  =  ^  =  {e)P  +  PA^(9)  +  (9)P  -  C.,{9)}  (D^{9)  +  (9))  (B  (9)P  -  C^(9))  +  el, 

(3.22) 

0=^  =  F^Q  +  QF^  -  rp(P-^)  +  ,  (3.23) 


where 


dP 


“1 


B^P  - 


Although  (3.21)-(3.23)  characterizes  extremals  to  (3.19),  we  have  not  yet  developed  a  reliable 
method  to  compute  these  extremals.  Note  that  standard  interior  point  descent  methods  (e.g., 
steepest  descent,  conjugate  gradient,  or  quasi-Newton  methods)  cannot  be  directly  applied  due  to 
the  nature  of  the  constraints.  For  example,  suppose  we  attempt  to  initialize  one  of  these  methods 
with  a  multiplier  (in  the  class  of  multipliers  for  the  given  uncertainty  set  )  represented  by  6o  and 
also  choose  an  initial  c  denoted  by  cq.  Then,  if  there  exists  a  positive-definite  solution  Pq  to  (3.15), 
the  REFP  is  solved  and  there  is  no  need  for  further  computations.  However,  suppose  there  is  no 
positive  definite  solution  Pq  to  (3.15).  Then,  Po)  cannot  be  used  to  initialize  an  interior  point 

descent  method  to  find  a  solution  to  the  optimization  problem  (3.19)  since  this  class  of  methods 
requires  an  initial  feasible  interior  point.  What  is  needed  is  a  numerical  technique  that  is  able  to 
find  a  solution  to  (3.b9)  by  starting  with  a  nonfeasible  point  (^o^Co,  Fo)*  This  is  accomplished  in 
the  next  section  using  a  probabililty-one  homotopy  algorithm. 


4.  Probability-One  Homotopy  Algorithms  for  Robust  Controller  Synthesis 

Homotopies  have  traditionally  been  studied  as  a  part  of  topology  and  have  found  significant 
application  in  nonlinear  functional  analysis  and  differential  geometry.  However,  it  has  only  been 
in  the  last  two  decades  that  homotopies  have  been  used  for  practical  numerical  computation.  The 
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algorithms  described  here  are  probability  one,  globally  convergent  homotopy  algorithms.^®’  This 
class  of  algorithms  is  related  to,  but  distinct  from,  other  homotopy  methods  such  as  continuation 
methods.^®  An  advantage  of  probability-one  algorithms  is  that  under  certain  conditions,  it  is  po-s- 
sible  to  guarantee  convergence  of  the  algorithm  from  an  arbitrary  starting  point  with  probability 
one.  (This  does  not  imply  that  the  algorithms  are  guaranteed  to  converge  to  global  optimum,  only 
that  they  can  converge  to  a  local  extremum  from  arbitrary  starting  points.)  These  algorithms  also 
relieve  some  of  the  technical  difficulties  sometimes  encountered  when  implementing  alternative  ho¬ 
motopy  methods.  For  example,  they  allow  the  zero  curve  of  the  homotopy  map  to  be  nonmonotonic 
in  the  homotopy  parameter  A,  and  avoid  singularities  along  the  zero  curve. 

Consider  a  function  F  :  xlZ-*  that  is  C^.  Given  7/  G  TJ,  it  is  desired  to  find  x  G 

such  that 

0  =  F(.t.7/).  (4.1) 


This  is  a  standard  zero  finding  problem.  In  the  context  of  the  robustness  analysis  results  of  the 
previous  section 


x  =  {0,e),  iV  =  9+l, 

^  fdc  dc\ 

W' 


(4.2) 

(4.3) 


and  corresponds  to  some  desired  lower  bound  on  the  multivariable  stability  margin.  Note  that 
0=1^  and  0  =  ^  are  implicitly  satisfied  by  choosing  P  as  the  (maximal)  solution  of  the  Riccati 
equation  (.3.22)  (or  (3.15))  and  Q  a,s  the  solution  of  the  Lyapunov  equation  (3.23). 


Let  .To  =  fo]  represent  a  feasible  multiplier,  a  stablizing  compensator  and  a  positive  c. 
Furthermore  let  70  be  chosen  small  enough  such  that  there  exists  a  positive-definite  solution  Pq  to 
(3.15).  (It  is  assumed  that  70  <  7/  such  that  the  robustness  problem  is  not  trivial.) 


We  let 


7(A)  =  (1  -  A)7o-f  A7/ 


and  define  the  probability-one  homotopy  map  p  :  [0, 1)  x  TZ^ 


by 


(4.4) 


/)(A,  t)  =  AF(t,  7(A))  +  (1  -  A)(t  -  .To).  (4.5) 

Obviously,  0  =  />  has  the  unique  solution  .tq  and  p  =  F(x.'fj).  These  are  necessary  conditions  for 
the  homotopy  map.  In  the  context  of  the  robustness  problem  of  Section  3,  this  homotopy  map  has 
the  desirable  property  that  it  can  be  initialized  with  any  feasible  multiplier.  In  addition,  for  any 
A  G  [0, 1)  the  corresponding  point  on  the  zero  curve  (t,  A)  corresponds  to  a  controller  and  multiplier 
that  guarantees  the  level  of  robustness  corresponding  to  7(A)  since  the  Riccati  equation  (3.15)  (or 
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(3.22))  with  the  constraints  (3.16)  are  satisfied  with  7  =  7(A).  Hence,  each  point  on  the  zero  curve 
(0  =  p(A.,r).  A  6  [0, 1)),  is  physically  meaningful  even  though  F(.t,7(A))  7^  0  for  0  <  A  <  1. 

4.1.  Probability-One  Homotopy  Algorithm 

Complete  details  of  the  numerical  algorithm  are  in  Reference  27.  A  sketch  follows. 

1.  Set  A  =  0,  .T  =  .To- 

2.  Evaluate  the  homotopy  map  p  and  the  .Jacobian  of  the  homotopy  map  Dp. 

3.  Predict  the  next  point  on  the  homotopy  zero  curve  using,  e.g.,  a  Hermite  cubic  inter- 
polant. 

4.  For  k  =  0, 1.2, ...  until  convergence  do 

^(fc+i)  g(k)  _  l£>p(Z(^^)fp{Z(% 

where  [Dp(Z)f  is  the  Moore-Penrose  pseudoinverse  of  Dp(Z).  Let  (ti,  Ai)  =  limjt_.co-^*^. 

•5.  If  Ai  <  1,  then  set  t  =  Ti,  A  =  Ai,  and  go  to  to  step  (2). 

6.  If  Ai  >  1,  compute  the  solution  x  at  A  =  1  using,  e.g.,  inverse  linear  interpolation.^' 

5.  The  Popov  Multiplier  and  an  Illustrative  Numerical  Example 

In  this  section  we  consider  the  special  case  of  the  Popov  multiplier.  In  particular  we  let  M($)  = 
Ns.,  where  H  and  N  are  real  and  diagonal.  This  multiplier  is  in  the  class  of  multipliers  for 
real  parametric  uncertainty  and  can  also  be  applied  to  complex  uncertainty  by  choosing  iV  =  0, 
in  which  case  H  is  the  constant  F-scale  matrix  considered  in  Reference  17,  38.  Note  that  if  M{s) 
is  the  Popov  multiplier  and  if  T!y(s)  is  asymptotically  stable,  then  M{s)T~f{s)  is  asymptotically 
stable.  Hence,  the  framework  of  Section  2  applies  with  Ma(s)  =  M{s)  and  Mb{s)  =  I.  Next  we 
use  the  Popov  multiplier  to  consider  the  case  of  nonrepeated,  real  parametric,  scalar  uncertainty 
(i.e.,  I,-  =  and  ki  =  1  for  i  =  1,  ...,p). 

5.1.  Nonrepeated,  Real  Parametric,  Diagonal  Uncertainty 

It  is  important  to  note  that  for  real  parametric  uncertainty,  D,  the  direct  feedthrough  term  in 
the  state  space  realization  of  G(s)  in  Figure  1,  is  always  zero.  Hence,  the  state  space  realization  of 


Ty(s)  (see  (2.13)  and  (2.14))  becomes 


Ty(s)  ~ 

A  +  jBC 

y/^B 

V^c 

I 

or.  equivalently. 

^s)  =  Tspls)  +  /, 

where  Tsp{s)  is  the  strictly  proper  transfer  function  given  by 

Tspis)  ~ 

A  +  'yBC 

v/275  ■ 

y/^C 

0 

(5.1) 


(5.2) 


(5.3) 


Note  that  for  nonrepeated,  real  parametric,  diagonal  uncertainty  the  Popov  multiplier  (with  H  and 
N  diagonal)  is  a  compatible  multiplier.  In  this  case  the  robustness  condition  (2,12)  becomes 

+  j(^N)(Tfip{jLo)  +  /)]  >0,  u  ^  TZU  oo, 


or,  equivalently, 

Ke[(H^  +  ju;N)Tsp(ju;)]  +  Ee[H^  +  ju;N]  >0,  u;  €  7^  U  oo. 
It  follows  from  Lemma  3  of  Reference  19  that 
iH^  +  Ns)T,pis) 


(5.4) 


(5.5) 


A  +  tBC 

y/¥fB 

y/^iH^C  +  NC(A  +  ^BC)) 

2'yNCB 

(5.6) 

Then,  since  HeffT^  +  jwiV]  =  =  He  the  robustness  condition  (5.4)  is  equivalent  to 

He  TwUlo)  >  0,  a;  €  7^  U  oo,  (5.7) 

where 


“  ■  ■  ~  s/^{H^C  +  iVC'(A+  tHC))  I  +  2'^NCB  J  ‘  ^  ' 

The  next  theorem  follows  immediately  from  (5.7),  (5.8),  and  Lemma  1. 

Theorem  6.  Let  J,  =  i  =  l,...,p,  assume  G{s)  is  asymptotically  stable,  and  define 

1.  If  there  exist  real  diagonal  H  and  N,P  >  0  and  c  >  0  such  that  NCB+B^C^  N+'y~^H^  >  0 
and 

0  =  A/P  +  PA^  +  €l 

+  [B^P  -  H^C  -  NCAyf[NCB  +  B^C'^N  +  j-^H^f^B'^P  -  H^C  -  NCA^l 

(5.9) 


or. 
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2.  if  there  exist  positive  definite,  diagonal  H^)  and  N .  P.  and  e  >  0  such  that 

-A^'^P  -  PA^.  -y/2^PB  +  +  ,/^AjC'^N  1 

-^B^P  +  y/^SC+y/^NCA^  2‘fNCBA2jNB'^C'^  +  2S-2€l  j  -  ’ 

(5.10) 

then  the  negative  feedback  interconnection  of  G'(.?)  and  A  is  asymptotically  stable  for  all 
A  e  A.,  . 


5.2.  Robust  Control  Synthesis  Using  the  Popov  Multiplier  for  a  Benchmark 
Problem 

To  illustrate  robust  control  synthesis  with  the  probability-one  homotopy  algorithm,  we  consider 
the  two-mass/spring  benchmark  system  shown  in  Figure  2  with  uncertain  stiffness  k.  A  control 
force  acts  on  the  body  1  and  the  position  of  body  2  is  measured,  resulting  in  a  noncolocated  control 
problem.  This  benchmark  problem  is  discussed  in  detail  in  Reference  39. 

The  open-loop  plant  (for  nii  =  =  1)  is  given  by 


Xp{t)  =  Ap(k)Xp{t)  +  Bpu{t)  +  Diw(t), 
y(t)  =  CpXp(t)  +  Z?2W.'(f), 
z(t)  =  EiXpit), 


where  2  =  .r2  is  the  performance  variable,  y  is  a  noise  corrupted  measurement  of  x^. 


Xp  =  [.ri ,  X2,  X3,  a;^] ,  Ap(k)  = 


0  0 
0  0 
0  0 
1  0 


Cp  =  Ei  = 


0 

0 

1 

0  ■ 

■  0  ■ 

0 

0 

0 

1 

.  Bp  = 

0 

-k 

k 

0 

0 

1 

k 

-k 

0 

0 

0 

[O  1 

0 

0 

o 

II 

‘1- 

We  desire  to  design  a  constant  gain  linear  feedback  compensator  of  the  form 


ic(<)  =  AcXc(t)  +  Bcy(t), 
u(t)  =  -CcXc{t), 


(.5.11) 

(5.12) 

(5.13) 


(5.14) 

(5.15) 


such  that  the  closed-loop  system  is  stable  for  0.5  <  k  <  2.0  and  for  a  unit  impulse  disturbance  at 
t  =  0,  the  performance  variable  z  has  a  settling  time  of  about  15  s  for  the  nominal  system  (with 
k  =  ^nom  ~  1)‘ 


We  begin  by  constructing  the  uncertainty  feedback  system  as  shown  in  Figure  1.  It  is  assumed 
that  k  =  k„om  +  AA:.  The  perturbation  in  Ap{k)  due  to  a  change  Ak  in  the  stiffness  element  k  from 
the  nominal  value  A’nom  is  given  by 


where  B,-,  = 


*dp(Ar)  -^plArnom)  —  AAp  —  B^AkCQ^ 
\T 


(5.16) 


0  0  1—1  ,  and  Co=Jl  —1  0  0  .  Assuming  negative  feedback,  the 


closed- loop  state  matrix  is  given  by 


A{k)  = 


Ap(k)  BpC'c 

BcCp  Ac 


(5.17) 


Next,  define  B  —  ^  Bq  04xi  ,  and  (7  =  Co  0ix4  so  that  G(s)  and  A  in  Figure  1  are  given 


by  G'(  .s) 


-^(^’nom) 

B  ' 

c 

0 

and  A  =  Ak. 


The  upper  bound  on  the  H2  cost  functional  discussed  in  Remark  7  is  given  by 


.7  =  ^tr[(F  -I-  2jC^NC)V], 


where  V  is  given  by 


V  = 


BcV^  BcV2BJ 


(5.18) 


(5.19) 


with  Vi  =  DiDi^ ,  V2  =  pD2D2^,  and  Vu  =  D\D2^ .  Also,  as  in  Remark  7,  the  Riccati  equation 

in  (5.9)  is  modified  slightly  in  that  c/  is  replaced  by  ei?  where  R  is  given  by 


R  = 


R\  —RuCc 

-CJR(2  CJR2Cc 


(5.20) 


where  Ri  =  R2  =  p,  and  R12  =  y/p  Ef.  Here  the  parameter  p  is  used  as  a  design  parameter 

to  increase  or  decrease  the  authority  of  the  controller. 

It  can  be  seen  that  the  diagonal  H  and  N  of  the  Popov  multiplier,  reduce  to  scalars  for  this 
particular  example.  The  variable  x  in  (4.2)  is  given  by 

\T 


X  = 


H  N  €  vec^(Ac)  vec^{Bc)  vec^(Cc)  ]  , 
and  consequently,  the  function  F  in  (4.3),  is  given  by 

^(®’T)=[fF  m  If  ^^€^(1^) 


(5.21) 


(5.22) 


The  initial  point  xq  is  chosen  in  the  following  manner.  Hq,  No,  and  cq,  are  chosen  arbitrarily 
as  10,  10,  and  1,  respectively.  The  initial  controller  {Acft,Bc,o,Cc,o)  is  an  LQG  controller  for  the 
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nominal  plant  corresponding  to  p  =  0.001.  No  robustness  is  expected  of  this  controller  and  hence 
7o  in  (4.4)  is  chosen  close  to  zero  (i.e.,  70  =  0.01). 

The  controller  transfer  function  obtained  by  the  probability-one  homotopy  algorithm  in  Sub¬ 
section  4.1  is  given  by 


H(s) 


2819  {s  +  0.2079  )[(.s  -  0.0978)^  +  0.806-3^] 

[is  4-  4.004)'^  +  1.82942][(s  +  3.4747)^  +  9.974.5'-*]  ’ 


(5.23) 


This  controller  is  guaranteed  by  the  theory  to  be  robust  for  the  range  0.5  <  k  <  2.0  and  this  was 
also  verified  by  a  direct  search.  The  settling  time  for  the  system  was  chosen  to  be  the  time  required 
for  the  displacement  of  mass  2  to  reach  and  stay  within  the  interval  [-0.1m,  0.1m].  The  controller 
is  seen  to  satisfy  the  settling  time  objective  when  connected  to  the  nominal  model  corresponding 
to  A’  =  1  N/m,  as  can  be  seen  from  the  impulse  response  of  the  closed- loop  system  in  Figure  3. 


6.  Conclusion 


This  paper  has  demonstrated  that  fixed-architecture,  robust  control  design  using  general  fixed- 
structure  multipliers  can  be  formulated  as  a  Riccati  equation  feasibility  problem,  a  nonlinear, 
algebraic  feasibility  problem.  Probability-one  homotopy  algorithms  were  proposed  to  solve  this 
feasibility  problem.  These  algorithms  differ  from  previously  developed  continuation  algorithms, 
developed  exclusively  for  the  case  of  the  Popov  multiplier,  in  that  they  can  be  initialized  with  any 
admissible  multiplier  and  stabilizing  compensator.  Like  other  probability-one  homotopy  algorithms 
they  also  tend  to  be  better  conditioned  than  the  alternative  continuation  algorithms.  The  results 
were  specialized  to  the  special  case  of  Popov  multipliers  and  fixed-architecture  robust  control  design 
was  illustrated  using  a  standard  benchmark  problem. 

The  key  step  in  the  practical  implementation  of  these  homotopy  algorithms  is  the  numerical 
implementation  of  the  computation  of  the  Jacobian  of  the  homotopy  map.  Analytical  expressions 
for  the  Jacobian  tend  to  be  complex  and  difficult  to  derive  due  to  this  complexity.  For  the  practical 
implementation  of  homotopy  algorithms  for  more  general  forms  of  the  multipliers  it  is  important 
to  develop  reliable  symbolic  matrix  differentiation  routines.  The  special  matrix  structures  that 
arise  in  the  computation  of  the  Jacobian  must  also  be  exploited  to  speed  up  the  algorithms.  An 
illustration  of  this  is  given  in  Reference  40. 
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Figure  1;  Standard  Uncertainty  Feedback  Configuration 
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Figure  2;  Two  Mass/Spring  System 
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Figure  3:  Impulse  Response  of  Closed-Loop  System 
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Abstract 

Reduced-order  control  is  important  in  control  engineering  practice  due  to  inevitable  limitations 
on  the  throughput  of  the  control  processors.  This  paper  considers  the  direct  design  of  H2  optimal, 
reduced-order  controllers  via  parameter  optimization  approaches.  Two  classes  of  numerical  methods 
are  considered:  1)  descent  methods,  and  2)  continuation  (or  homotopy)  methods.  In  particular, 
four  standard  descent  methods  are  compared  with  a  recently  developed  continuation  algorithm  by 
using  three  examples  appearing  in  the  literature.  The  continuation  algorithm  is  seen  to  be  much 
more  reliable  than  the  descent  methods. 

1.  Introduction 

One  of  the  deficiencies  of  modern  control  laws,  developed  by  simply  solving  a  pair  of  decoupled 
Riccati  equations,  in  particular,  linear-quadratic-guassian  (LQG)  control  and  full-order  Hoo  control, 
is  that  the  resultant  control  laws  are  always  of  the  order  of  the  design  plant.  These  techniques, 
though  relatively  easy  to  implement  computationally,  do  not  allow  the  designer  to  constrain  the 
architecture  (e.g.  order  and  degree  of  centralization)  of  the  controller.  Such  constraints  are  often 
necessary  in  engineering  practice  due  to  throughput  limitations  of  the  control  processors.  Reduced- 
order  control  is  therefore  of  paramount  importance  in  practical  control  design.  Tliis  paper  focuses 
on  the  design  of  H2  optimal,  reduced-order  controllers. 

Two  main  approaches  have  been  developed  to  solve  the  H2  optimal,  reduced-order  design  prob¬ 
lem.  The  first  approach  attempts  to  develop  approximations  to  the  optimal  reduced-order  controller 


by  reducing  the  dimension  of  an  LQG  controller(Yousuff  and  Skelton  1984a,  YousufF  and  Skelton 
1984b,  Anderson  and  Liu  1989,  Villemagne  and  Skelton  1988,  Liu  et  al.  1990).  These  methods 
are  attractive  because  they  require  relatively  little  computation  and  should  be  used  if  possible. 
Unfortunately,  they  tend  to  yield  controllers  that  either  destabilize  the  system  or  have  poor  perfor¬ 
mance  as  the  requested  controller  dimension  is  decreased  or  the  requested  control  aufliority  level 
is  increased.  Hence,  if  used  in  isolation,  these  methods  do  not  yield  a  reliable  methodology  for 
reduced-order  design.  In  addition,  these  methods  do  not  extend  to  the  design  of  decentralized 
controllers.  However,  it  should  be  mentioned  that,  in  regards  to  reduced-order  control  design, 
the  indirect  approaches  at  worst  are  valuable  in  providing  good  initial  conditions  for  the  direct 
approaches  described  below. 

In  contrast  to  controller  reduction,  direct  approaches  attempt  to  directly  synthesize  an  optimal, 
reduced-order  (or  decentralized)  controller  by  a  numerical  optimization  scheme.  There  are  two  main 
classes  of  parameter  optimization  approaches  to  direct  control  design.  The  first  class  relies  on  the 
use  of  descent  methods  (Kramer  and  Calise  1987,  Kuhn  and  Schmidt  1987,  Kwakernaak  and  Sivan 
1972,  Ly  et  al.  1985,  Mukhopadhyay  1982,  Mukhopadhyay  1987,  Voth  and  Ly  1991).  Algorithms 
in  tliis  class  reduce  the  cost  at  each  iteration.  The  second  class  reUes  on  the  use  of  continuation 
methods  (Collins  et  al.  1994,  Mercadal  1991).  In  contrast  to  the  descent  methods,  the  H2  cost 
is  not  necessarily  reduced  at  each  iteration.  It  should  be  mentioned  that  continuation  algorithms 
(Collins  et  al  1996b)  have  also  been  developed  to  solve  the  “optimal  projection  equations,”  a  set 
of  four  coupled  Lyapunov  and  Riccati  equations  that  characterize  the  H2  optimal,  reduced-order 
compensator.  However,  this  approach  will  not  be  considered  here. 

From  a  practical  design  perspective  it  is  important  to  determine  which  class  of  methods  tends 
to  be  more  numerically  robust.  As  with  the  vast  majority  of  numerical  methods  tor  nonconvex  opti¬ 
mization  problems,  answers  to  these  questions  are  extremely  difiicult  to  prove  analytically.  Instead, 
we  must  rely  on  numerical  experimentation  to  observe  trends.  Hence,  in  this  paper  the  behavior  of 
some  standard  descent  methods  (i.e.,  steepest  descent,  conjugate  gradient,  BFGS  Quasi-Newton, 
and  Newton’s  method)  (Fletcher  1987)  are  compared  to  the  corresponding  behavior  of  the  continu¬ 
ation  algorithm  of  (Collins  et  al.  1994)  by  considering  design  for  three  reduced-order  control  design 
problems  appearing  in  the  literature.  The  results  clearly  indicate  that  the  continuation  algorithm 
tends  to  be  more  numerically  robust  and  is  most  efficient  when  the  controller  is  constrained  to  a 
tridiagonal  form. 

The  paper  is  organized  as  follows.  Section  2  formulates  the  H2  optimal,  reduced-order  dynamic 
compensation  problem  as  a  constrained  parameter  optimization  problem  and  discusses  various  basis 
options  for  the  controller  in  order  to  reduce  the  size  of  the  controller  parameter  vector.  Section 


3  briefly  describes  the  algorithms  for  the  descent  and  continuation  methods  for  the  H2  optimal, 
reduced'Order  control  problem.  The  descriptions  emphasize  how  the  constraint  that  the  controller 
be  stabilizing  is  taken  into  account  in  the  algorithms.  Section  4  uses  the  three  examples  mentioned 
above  to  present  a  comparison  of  both  the  numerical  robustness  and  speed  of  convergence  of  the 
descent  and  continuation  methods.  Finally,  Section  5  presents  conclusions. 

Notation 


n  X  1  real  vector 
transpose  of  z 

tr(M )  trace  of  the  square  matrix  M 

E  expectation  operator 

dom{g)  domain  of  the  function  g(-) 


2.  H2  Optimal,  Reduced-Order  Dynamic  Compensation 
2.1.  Problem  Formulation 
Consider  the  system 

x{t)  =  Ax{t)  +  Bu{t)  +  Diw{t)  (2.1) 

y{t)  =  Cx{t)  +  Du{t)  +  D2w{t)  (2.2) 

z{t)  =  Eix{t)  +  E2u{t)  (2.3) 

where  x  €  7^"*,  «  €  7^"“,  y  €  ^  €  7^"*,  w  €  7^”“  is  white  noise  with  unit  intensity,  D2  has 

full  row  rank,  and  E2  has  full  column  rank.  We  desire  to  design  a  order  dynamic  compensator, 

Xc(t)  =  AcXc(t)  +  Bcy{t)  (2.4) 

u{t)  =  -CcXc{t)  (2.5) 

where  ric  <  which  minimizes  the  steady  state  performance  criterion 

J{AcyBc,Cc)  =  ^  E{z^{t)z{t)\.  (2.6) 

In  the  frequency  domain,  the  cost  in  (2.6)  may  be  interpreted  as 

1 

J(A„  Be,  Cc)  =  ^J  ’  (2-7) 

where  G2«;(s)  is  the  transfer  function  from  to  to  z  and  the  right  hand  side  is  the  square  of  the  H2 
norm  of  Gzw{s)- 


The  state-space  evolution  of  the  closed-loop  system  corresponding  to  (2.1)-(2.5)  is  described  by 


where 


x(0  = 

Now,  the  cost  (2.6)  can  be  expressed  as 


a;(t) 

i  = 

®c(0 

i(/.)  =  Ax{t)  -f  Dw{i) 


A  -BCc 
BcC  Ac-BcDCc 


,  D  = 


Di 

BcD2 


J(Ac,5c,Cc)  =  ^Uin 


(2.8) 

(2.9) 

(2.10) 


where 


Ri  RuCc 


,  Ri^Ex^Ei,  Ru  =  2Ei^E2,  R2  =  E2^E2.  (2.11) 


C'Wa  CJR2Cc 

(Note  that  since  E2  has  full  column  rank,  R2  >  0.) 

To  guarantee  that  the  cost  J  is  finite  and  independent  of  initial  conditions  we  restrict  our 
attention  to  the  set  of  stabilizing  compensators,  Sc  =  {(AciBcCc)  :  A  is  asymptotically  stable}. 
Assume  {Ac,Bc,Cc)  €  Sc  and  define  Q  6  to  be  the  closed-loop  steady-state 

covariance,  i.e., 

0:=AQ  +  QA^  +  V  (2.12) 


where 


F  = 


Vi2B^ 
BcVT  BcV2B'^ 


,  Fi2  =  2T>i^£»2,  V2^D2^D2.  (2.13) 

(Note  that  since  D2  has  full  row  rank,  V2  >  0.)  The  cost  function  J  car.  now  be  expressed  as 

J{Ac,Bc,Cc,Q)^iTQR.  (2.14) 

The  objective  is  to  minimize  the  cost  function  J  subject  to  the  constraint  (2.12). 

The  Lagrangian  C  is  defined  by 

£(Ac,  Be,  Cc,  Q,  P)  =  tiQR  +  tT[PiAQ  -t-  QA^  +  V)]  (2.15) 

where  P  is  the  Lagrange  multiplier  matrix.  The  compensator  (Ac,  Bc^  Cc)  is  optimal  if  it  satisfies 
the  stationary  conditions 

^  >1  # '  w  /  - 

(2.16) 


J^_0  J^-0  ^  =  0 

dAc~  '  dBc~  '  dCc  ’ 


dQ 


=  Ap  +  pa^  +  r  =  o, 


(2.17) 


3 


dP 

Both  the  descent  and  continuation  algorithms  aim  at  finding  (Ac,Bc,Cc)  €  Sc  that  satisfy  the 
above  conditions. 

Subsequently,  we  will  represent  the  controller  by  a  parameter  vector  0,  for  example, 


0  = 


vec(Ac) 

vec(Bc) 

vec(Cc) 


(2.19) 


Let  the  mapping  from  a  state  space  representation  of  a  controller  (Ac,  Bci  Cc)  to  the  parameter 
vector  0  be  given  by  g{-),  such  that 

0  =  g{Ac,Bc,Cc)  (2.20) 

and  define 

0  =  =  g{Ac,  Be,  Cc) :  (Ac,  B^  Cc)  €  Sc  U  dom(5)}  (2.21) 

Now,  assuming  ^  €  0,  the  H2  cost  functional  and  the  corresponding  Lagrangian  can  be  expressed 
respectively  as  J(0,Q)  and  C{0,Q,P).  The  problem  is  therefore  to  find  0  €.  Q  such  that 

O=^i0,Q,P).  (2.22) 


subject  to  (2.18)  and  (2.17), 


2.2*  Reduction  of  the  Dimension  of  the  Controller  Parameter  Vector  (0) 

It  is  desired  that  the  parameter  vector  6  be  as  small  as  possible.  Hence,  we  desire  to  represent 
the  controller  matrix  with  the  fewest  parameters  possible.  The  minimal  number  of  parameters  Pmin 
with  which  a  compensator  can  be  represented  is  given  by  (Denery  et  al.  1971,  Martin  and  Bryson 
1980) 

Pimn=nc{m  +  l).  (2.23) 

One  canonical  form  which  allows  representation  of  a  controller  with  a  minimal  number  of  pa¬ 
rameters  is  the  modal  form  described  in  (Martin  and  Bryson  1980).  This  form  will  be  called  here 
the  Second-Order  Polynomial  (SP)  form.  For  this  parameterization  a  triple  {Ac,Bc,Cc)  of  order 
He  has  the  following  structure.  ,  , 


Ac  =  block-diag{Ac,i,  Ac,2,...,Ac,nJ 


(2.24) 


where  Ac,i  is  2  X  2,  t  =  {1, 2, 7ic}  arid  each  Ac,i  (with  the  exception  of  if  the  row  dimension 
of  Ac  is  odd)  has  the  form 


•'4c,t  — 


0  1 
4?  j 


(2.25) 


to  allow  for  either  a  complex  conjugate  set  of  poles  or  two  real  poles.  Be  is  completely- full  and 

Cc  =  lCc,i,Cc,2,,-,Cc,rh  (2-26) 


where  (7c,t  the  form 


1  0 
*  * 

*  * 


(2.27) 


The  controller  canonical  form  described  in  (Kailath  1980)  also  allows  representation  of  a  con¬ 
troller  with  a  minimal  number  of  parameters.  For  SISO  systems  in  controller  canonical  form  the 
Ac  matrix  is  a  companion  matrix.  In  particular,  Ac  has  the  form 


Ac^ 


0  10  0 

0  0  10 

0  0  0  1 

*  +  +  + 


0 

0 

0 

0 

* 


(2.28) 


In  addition, 


Bc  = 


(2.29) 


and  Cc  is  completely  full.  A  dual  form  of  the  controller  canonical  form  is  the  observable  canonical 
form  (Kailath  1980). 

Another  minimal  parameter  basis  is  the  ‘‘input  normal  Riccati  form”  of  (Davis  et  cL  1994).  In 
this  basis  Ac  is  completely  determined  by  Be  and  Cc- 

It  is  also  possible  to  represent  the  controller  in  a  basis  where  the  number  of  free  parameters  p 
satisfies 

Pmin  <  P<  Pmax  =  Re  (nc  +  m -h /).  -  (2.30) 

One  such  basis  is  the  tridiagonal  basis  (Geist  1991,  Parlett  1992)  in  which  the  controUer  state 
matrix  is  constrained  to  have  nonzero  elements  only  on  the  diagonal,  the  super-diagonal,  and  the 


sub-diagoiial.  That  is, 


+  * 

*  *  *  0 

,  (2-31) 

0  ■■■ 

*  *  _ 

and  Be  and  Cc  are  completely  fuU.  For  this  form  the  number  of  free  parameters  is  given  by 

P  —  Pram  d* 

It  is  important  to  recognize  that,  given  a  particular  basis,  there  are  stabilizing  controllers  that 
have  no  state  space  reaUzation  in  that  basis,  such  that  if  <7(-)  in  (2.20)  requires  the  representation 
of  a  controUer  in  a  given  basis  Sc  U  dom(5)  C  Sc-  Hence,  the  set  0,  defined  by  (2.21)  corresponds 
to  a  smaller  set  of  controllers  than  the  set  Sc-  When  the  controller  is  restricted  to  a  particular 
basis  in  the  algorithms  to  follow,  this  reduction  in  the  size  of  the  feasible  set  sometimes  leads  to 
numerical  ill-conditioning  or  even  algorithm  failure. 


3.  Parameter  Optimization  Algorithms 

This  section  first  gives  a  general  description  of  the  algorithms  corresponding  to  the  descent 
methods.  It  then  briefly  describes  a  continuation  algorithm.  Particular  attention  is  given  to  the 
modification  of  these  algorithms  to  take  into  account  the  constraint  0  6  0. 

3.1.  Descent  Methods 

Descent  methods  are  designed  to  search  for  solutions  to  the  unconstrained  optimization  problem 

imnJ(^).  (3-1) 

The  user  is  required  to  supply  an  initial  parameter  vector  Oq.  A  descent  algorithm  then  has  the 
following  structure. 

A  Descent  Algorithm 

1.  Let  k  =  0.  . 

2.  Determine  a  search  direction  dk- 

3.  Use  a  one  dimensional  line  search  to  find  a*  that  minimizes  J{6k  +  adk)  with  respect  to  a. 


4.  Set  Ok+i  =  Ok  +  oikdk 


5. 


If  the  gradient  ^{Ok+i)  is  sufficiently  small,  then  let  the  optimal  solution  0*  =  0k+i  and  stop, 


else  let  A;  =  fc  +  1  and  go  to  Step  2. 


Alternative  descent  methods  differ  primarily  in  the  way  they  compute  the  descent  direction  dk~ 


For  example,  in  the  steepest  descent  method  dk  corresponds  to  the  negative  of  the  gradient.  Con¬ 
jugate  gradient  and  Quasi-Newton  methods  compute  dk  using  only  cost  and  gradient  information 
while  Newton’s  method  requires  computation  of  the  Hessian  matrix.  Note  that  for  the  H2  optimal. 


reduced-order  control  problem  it  is  not  difficult  to  show  that  if  (2.18)  and  (2.17)  are  satisfied,  then 


the  graxlient  satisfies 


^_dC 

'm~  d0' 


(3.2) 


Hence,  the  gradient  may  be  computed  by  constructing  and  differentiating  the  Lagrangian. 

Recognize  that  the  H2  optimal,  reduced-order  control  problem  is  not  the  unconstrained  opti¬ 
mization  problem  (3.1)  but  is  actually  the  constrained  optimization  problem 


minJ(0)  (3-3) 

eg© 

where  0  is  defined  by  (2.21).  One  way  to  take  into  account  the  constraint  0  G  0  is  to  modify  the 
line  search  subalgorithm  of  Step  3  to  ensure  that  if  G  0,  0k+i  is  also  in  0.  (It  is  assumed  that 
00  G  0.)  An  example  of  such  a  modification  is  discussed  in  (Kuhn  and  Shmidt  1987).  The  descent 
algorithms  compared  in  this  paper  use  the  following  modification  to  the  line  search  algorithm  of 
(Fletcher  1987).  During  the  bracketing  phase  of  the  line  search  algorithm  a  stability  check  is  carried 
out  to  make  .sure  that  the  right  end  of  the  bracket  results  in  a  stable  closed-loop  system.  If  not, 
then  the  bracket  is  halved.  This  simple  modification  ensures  that  aU  points  within  the  bracket 
result  in  a  stable  closed-loop  system. 


®  3.2.  Continuation  Methods 

Continuation  techniques  can  be  used  to  solve  the  zero  finding  problem 

0  =  f{0),  (3.4) 

where  /  :  TV  TV.  In  the  context  of  optimal,  reduced-order  control,  (3.4)  corresponds  to 
(2.22).  Continuation  techniques  work  by  finding  a  function  H  :  TV  X  [0, 1)  TV  that  satisfies 
certain  properties,  including  the  following: 


7 


2.  0  =  H(0,0)  has  an  easily  found  or  known  solution  ^o- 


They  then  trace  the  zero  curve  described  by 

0  =  H{e,X),  A€[0,1).  .  (3-5) 

Tliis  is  accomplished  by  differentiating  (3.5)  with  respect  to  A  to  obtain  Davidenko’s  differential 
equation 

0  =  Hx{9,X)  +  He{e,X)9x{X)  (3-6) 

where  Hx  =  Hg  =  ,  and  9x  =  ,  which  together  with  ^(0)  =  9o  defines  an  initial  value 

problem.  Predictor-corrector,  numerical  integration  schemes  are  then  used  to  solve  this  initial  value 
problem,  that  is  to  follow  the  curve  (3.6)  from  the  solution  of  0  =  H{9,0)  to  a  solution  9*  of 
0  =  1).  In  particular  a  continuation  algorithm  has  the  following  structure. 

A  Continuation  Algorithm 


1.  Let  A  =  0  and  ^(A)  =  9o. 

2.  Use  (3.6)  to  compute  the  tangent  vector  9x,  such  that  ^a(A)  =  -Hg{9,X)~^JIx{9,X). 

3.  For  some  AA  such  that  A  =  A  +  AA  <  1,  use  current  and  past  values  of  H  and  Hx  to  predict 
9(X  +  AA)  by  using  polynomial  curve  fitting. 

4.  Let  A  <—  A  +  AA  and  9o  be  the  prediction  of  0(A). 

5.  For  k  =  0, 1,2, ...  until  convergence,  do 

^k+i  =6k-  Hg{9k,  X)~^9k. 

Then,  let  0(A)  =  0fc+i. 

6.  If  A  <  1,  go  to  Step  2,  else  if  A  =  1,  then  let  the  solution  0*  =  0(A)  and  stop. 

The  initializing  controller  0o  in  the  algorithm  for  H2  optimal,  reduced-order  control  (Collins 
et  al  1994,  Collins  et  al.  1996a)  is  usually  found  by  applying  a  controller  reduction  method 
such  as  balanced  controller  reduction  (Yousuff  and  Skelton  1984a)  to  ^  low  authqrity  LQG  con¬ 
troller  since  this  usually  yields  a  nearly  optimal,  reduced-order  controller.  The  initial  weights 
(£1)0,  (ili2)o,  (1^2)0,  (Vi)o,  (Vi2)o,  (V2)o  corresponding  to  the  low  authority  LQG  controUer  are 


then  deformed  into  the  desired  weights  along  the  homotopy  path.  The  reader  is  referred  to  (Collins 
et  al.  1994)  for  further  details. 

The  algorithm  of  (Collins  et  al.  1994)  also  assumes  that  the  prediction  ^(A  +  AA)  €  0  such 
that  it  corresponds  to  a  controller  that  stabilizes  the  closed-loop  system.  If  ^(A  H-  AA)  ^  0,  then 
the  algorithm  reduces  the  size  of  AA.  In  particular,  AA  *—  ^AA.  This  is  a  small  but  necessary 
modification.  As  subsequently  discussed,  it  is  at  this  point  that  the  algorithm  sometimes  falls. 


4.  Numerical  Examples 

4.1.  Description  of  Problems 

The  first  problem  is  a  noncoUocated  axial  vibration  control  problem  involving  an  axial  beam 
with  four  circular  disks  attached.  This  problem  was  introduced  in  (Cannon  and  Rosenthal  1984) 
and  also  studied  in  (Collins  et  al  1994).  The  plant  is  8th  order  while  we  consider  the  design  of  a 
4th  order  controller. 

The  second  problem  was  introduced  in  (Ly  et  al.  1985)  and  involves  flight  control  for  a  NAVION 
aircraft.  The  model  is  7th  order  and  we  consider  the  design  of  a  4th  order  controller. 

The  third  problem  was  introduced  in  (Martin  and  Bryson  1980)  and  involves  vibration  control 
of  a  flexible  spacecraft.  The  model  is  6th  order  while  we  again  consider  the  design  of  a  4th  order 
controller. 

The  system  matrixes  (A,  B,  C,  D)  and  the  weighting  matrices  (Ri,  R^,  Ru,  Vi,  V2,  V12)  are  given 
in  the  Appendix.  Note  that  the  matrices  R2  and  V2  are  multiplied  by  p  which  is  allowed  to  change 
from  10  to  1  in  order  to  deform  the  low  authurity  controller  to  a  higher  authority  controller  using 
the  homotopy  algorithm.  In  the  case  of  the  descent  algorithms  p  is  fixed  at  1. 

For  each  example,  a  low  authority  optimal  LQG  controller  (corresponding  to  p  =  10)  was  first 
designed.  The  order  of  this  controller  was  then  reduced  using  the  modified  balanced  controller 
reduction  technique  of  (Yousuff  and  Skelton  1984a).  This  reduced  order  sub-optimal  controller  was 
then  converted  into  an  optimal  low  authority  controller  using  a  few  Newton  iterations.  This  con¬ 
troller  was  used  as  the  starting  point  for  both  the  continuation  and  descent  optimization  methods. 
These  numerical  examples  are  included  in  Tables  1  through  3  and  were  run  on  a  90  MHz,  Pentium 
PC.  Table  4,  which  only  presents  data  for  the  continuation  algorithm,  was  produced  using  a  120 
MHz  Pentium  PC. 

We  also  design  higher  order  controllers  for  all  three  examples  using  continuation  algorithms 
with  the  controller  unconstrained  and  with  the  controller  constrained  to  the  tridiagonal  basis,  in 


order  to  compare  the  two  methods.  These  numerical  examples  are  included  in  Table  4  and  weie 
run  on  a  120  MHz  Pentium  PC. 


Method /Basis 

Unconstrained 

Tridiagonal 

SPF 

CCl 

M  Flops 

Time 

(sec) 

M  Flops 

Time 

(sec) 

M  Flops 

Time 

(sec) 

M  Flops 

Time 
'  (sec) 

Continuatioii 

19.8 

55 

14.8 

50.8 

374 

1422 

349.8 

1987 

Newton 

13.85 

41.1 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

BFGS 

13.45 

41.3 

3.4 

15.3 

f-ls 

f-ls 

f-ls 

f-ls 

Conjugate  Grad. 

56.5 

167.5 

13.7 

46.8 

f-osc 

f-osc 

f-ls 

f-ls 

Steepest  Des. 

227 

727 

68.3 

207.3 

f-osc 

f-osc 

f-osc 

f-osc 

Table  1:  FOUR  DISKS 


4*2.  Observations 

Comparison  of  the  various  algorithms  are  given  in  Tables  1,  2,  and  3.  The  letter  ‘T”  denotes 
failure  of  an  algorithm.  In  particular  ‘T-ls”  denotes  failure  due  to  the  step  length  parameter 
computed  by  the  line  search,  becoming  extremely  small.  This  occurs  because  of  the  modification  to 
the  line  search  algorithm,  to  make  sure  that  the  controller  results  in  a  stable  closed-loop.  During 
the  line  search  if  the  resulting  controller  violates  this  constraint,  as  previously  described,  the  step 
size  is  reduced  by  half.  This  sometimes  results  in  extremely  small  step  sizes  for  which  there  is 
no  appreciable  decrease  in  cost  or  change  in  the  search  direction.  Hence  the  algorithm  effectively 
comes  to  a  stop. 

Similarly,  "f-A”  denotes  failure  due  to  the  nei  d  for  extremely  small  increments  in  the  homotopy 
parameter  A.  This  occurs  because  of  the  modification  to  the  prediction  step  in  the  continuation 
algorithms  to  avoid  unstable  closed-loops.  If  during  any  of  these  steps,  an  unstable  closed  loop  is 
obtained,  then,  as  previously  described,  the  increment  to  the  continuation  parameter  A,  is  decreased 
by  half.  This  sometimes  results  in  the  need  for  extremely  small  increments  in  the  homotopy 
parameter  A  for  which  there  is  no  appreciable  change  in  the  controller  parameters 

The  indicator  “f-osc”  denotes  failure  in  a  descent  algorithm,  due  to  oscillatory  behavior  of 
the  gradient  of  the  cost  function  as  the  algorithm  progresses.  Several  hundreds  of  iterations  are 
performed  without  any  perceptible  change  in  the  cost  function.  This  phenomena  was  observed 
primarily  in  the  steepest  descent  algorithm  and  is  a  well  known  deficiency  of  this  method  (Fletcher 
1987). 

Tables  1  through  3  reveal  that  constraining  the  basis  of  the  controller  led  to  increased  failure  in 


Metliod/Basis 

SPF 

CCl 

7 

Time 

(sec) 

M  Flops 

Time 

(sec) 

Continuation 

537 

HO 

BEBI 

f-A 

hb 

f-A 

f-A 

Newton 

f-ls 

mm 

rnEm 

HO 

f-ls 

f-ls 

f-ls 

f-ls 

BFGS 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

Conjugate  Grad. 

f-ls 

■iEM 

f-ls 

f-ls 

HSHI 

f-ls 

Steepest  Des. 

f-osc 

f-osc 

f-osc 

mam 

msam 

f-ls 

Table  2:  NAVION 


Method/Basis 

Unconstrained 

Tridiagonal 

SPF 

CCF 

M  Flops 

Time 

(sec) 

M  Flops 

Time 

(sec) 

M  Flops 

Time 

(sec) 

M  Flops 

Time 

(sec) 

Continuation 

10.38 

21.7 

21.1 

25.1 

23.4 

121.8 

94.85 

Newton 

11.7 

5.9 

9.88 

5.22 

7.25 

7.57 

8.9 

BFGS 

f-ls 

■H 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

f-ls 

Conjugate  Grad. 

25.3 

253.4 

159.8 

51.19 

517.7 

302.6 

Steepest  Des. 

HIOI 

147.2 

f-osc 

f-osc 

6.92 

285 

193 

Table  3;  SPACECRAFT 


Example(Controller  Order)/Basis 

Unconstrained 

Tridiagonal 

M  Flops 

Time 

(sec) 

M  Flops 

Time 

(sec) 

Fourdisks  (6‘*^  order) 

1028.5 

294.5 

412.8 

164 

Fourdisks  (8*^**  order) 

596 

809.9 

239.4 

Navion  (6**^  order) 

600.8 

1566.2 

457.7 

Spacecraft  (6**^  order) 

347.9 

130.1 

196.5 

98.4 

Table  4;  Continuation:  Unconstrained  vs  Tridiagonal 


both  the  continuation  and  descent  algorithms  (i.e.  by  both  the  “f-A”  and  the  f-ls  failure  modes). 
One  apparent  reason  for  this  is  that  constraining  the  controller  to  a  particular  basis  reduces  the 
feasible  subspace,  as  discussed  at  the  end  of  Section  2.  Algorithm  failure  due  to  basis  constraints 
was  also  observed  in  (Kuhn  and  Schmidt  1987). 

The  numerical  conditioning  of  the  algorithms  when  using  the  tridiagonal  basis  was  better  than 
when  using  the  second  order  polynomial  form  (SPF)  and  the  controller  canonical  form  (CCF)  and 
is  apparently  due  to  the  fact  that  the  tridiagonal  form  is  a  more  general  representation  than  SPF 
and  CCF.  In  fact,  SPF  is  a  special  case  of  a  tridiagonal  form.  This  phenomena  can  be  seen  by 
a  comparison  of  execution  times  for  the  respective  bases  in  Tables  1  through  3.  For  most  cases, 
the  execution  times  for  the  SPF  and  the  CCF  bases,  are  several  times  larger  than  the  those  for 
the  tridiagonal  basis,  even  though  the  SPF  and  CCF  are  minimal  parameter  bases,  while  the 
tridiagonal  is  not,  and  hence  have  smaller  parameter  vectors  0.  This  is  due  to  the  decreased 
numerical  conditioning  associated  with  the  use  of  the  SPF  and  CCF  basis. 

The  conjugate  gradient  method,  as  might  be  expected,  usually  converged  faster  than  the  steepest 
descent  method.  The  convergence  times  for  the  Newton,  BFGS  and  continuation  methods  were 
equivalent  and  usually  much  less  than  the  convergence  times  for  the  conjugate  gradient  and  the 
steepest  descent  methods.  However  aU  the  descent  methods  including  Newton  and  BFGS  methods 
fail  much  more  often  than  did  the  continuation  methods.  This  is  especially  true  when  the  controller 
basis  is  unconstrained  or  tridiagonal,  in  which  case  the  continuation  method  never  failed. 

In  Tables  1  through  3,  it  is  seen  that  the  run  times  of  the  continuation  algorithm  for  the  uncon¬ 
strained  and  tridiagonal  cases  were  very  similar.  However,  as  the  controller  dimension  increases, 
the  size  of  the  parameter  vector  associated  with  the  unconstrained  “basis”  increases  much  more 
rapidly  than  the  parameter  vector  associated  with  the  tridiagonal  basis.  Hence,  it  is  expected  that 
the  convergence  times  for  the  tridiagonal  case  will  increase  much  less  rapidly  than  the  unconstrained 
case  as  the  controller  dimension  increases.  This  is  confirmed  in  Table  4  which  was  generated  using 
a  120  MHz  Pentium  PC. 

5.  Conclusions 

This  paper  has  used  three  examples  to  compare  the  behavior  of  four  standard  descent  algo¬ 
rithms  with  a  recently  developed  continuation  algorithm  for  H2  optimal,  reduced-order  design. 
The  results  clearly  indicate  that  the  continuation  algorithm  is  much  more  numerically  robust  than 
any  of  the  descent  algorithms,  especially  when  the  controller  basis  is  unconstrained  or  tridiagonal. 
The  numerical  conditioning  of  the  algorithms  always  decreased  when  the  controller  was  constrained 


to  either  the  second  order  polynomial  form  or  the  controller  canonical  form  which  are  minimal 
parameter  bases.  However  numerical  conditioning  of  the  algorithm  when  using  the  tridiagonal 
form,  a  nonminimal  parameter  basis,  was  better  than  when  using  a  minimal  parameter  basis. 
The  numerical  conditioning  of  the  continuation  algorithms,  when  using  the  tridiagonal  basis  was 
unchanged  or  at  worst  marginally  decreased,  as  compared  to  the  unconstrained  case.  Hence,  when 
using  the  tridiagonal  basis,  the  advantage  of  a  smaller  parameter  vector  usually  outweighed 
the  disadvantage  of  reduced  numerical  conditioning  due  to  a  basis  constraint.  This  advantage,  as 
expected,  became  more  apparent  as  the  order  of  the  controller  was  increased. 
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APPENDIX 


The  state  space  descriptions  of  the  3  examples  are  as  follows. 
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Cost-Effective  Parallel  Processing  for 
H^/  Controller  Synthesis 

Yuzhen  Ge+  Layne  T.  Watson^  Emmanuel  G.  Collins,  Jr.* 

Abstract 

A  distributed  version  of  a  homotopy  algorithm  for  solving  the  H^/H^  mixed- 
norm  controller  synthesis  problem  is  presented.  The  main  purpose  of  the  study  is  to 
explore  the  possibility  of  achieving  high  performance  with  low  cost.  Existing  UNIX 
workstations  running  PVM  (Parallel  Virtual  Machine)  are  utilized.  Only  the  Jaco¬ 
bian  matrix  computation  is  distributed  and  therefore  the  modification  to  the  original 
sequential  code  is  minimal.  The  same  algorithm  has  also  been  implemented  on  an 
Intel  Paragon  parallel  machine.  Our  implementation  shows  that  acceptable  speedup  is 
achieved  and  the  larger  the  problem  sizes,  the  higher  the  speedup.  Comparing  with  the 
results  from  the  Intel  Paragon,  the  study  concludes  that  utilizing  the  existing  UNIX 
workstations  can  be  a  very  cost-eflfective  approach  to  shorten  computation  time.  Fur¬ 
thermore,  this  economical  way  to  achieve  high  performance  computation  can  easily  be 
realized  and  incorporated  in  a  practical  industrial  design  environment. 

1  Introduction 

H^/H^  mixed- norm  controller  synthesis  is  an  important  and  interesting  technique  in  mod¬ 
ern  control  design  which  provides  the  means  for  simultaneously  addressing  and 
performance  objectives.  In  practice  such  controllers  provide  both  nominal  performance 
(via  suboptimal  H^)  and  robust  stability  (via  Hence  mixed-norm  synthesis  provides 

a  technique  for  trading  off  performance  and  robustness,  a  fundamental  objective  in  control 
design. 

The  mixed- norm  problem  has  been  addressed  in  a  variety  of  settings.  One 

treatment  utilized  an  cost  bound  as  the  basis  for  an  auxiliary  nonconvex  constrained 
minimization  problem,  which  is  very  difficult  to  solve  without  the  global  convergence  of  ho¬ 
motopy  methods.  A  successful  homotopy  algorithm  based  on  the  Ly  form  parametrization 
has  been  developed  (Ge  et  al  1994). 

The  control  design  algorithms  will  be  used  for  controller  design  of  systems 

such  as  the  four  disk  system  of  (Cannon  and  Rosenthal  1984).  This  system  is  especially 
representative  of  the  type  of  vibration  control  problems  that  arise  in  industrial  problems 
involving  rotating  turbomachinery.  H^/H^  design  will  be  used  to  develop  controllers  that 
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are  robust  with  respect  to  unmodeled  dynamics  and  also  guarantee  a  certain  measure  of 
nominal  performance. 

It  should  be  mentioned  that  theory  has  been  used  in  (Haddad  et  al.  1993)- 

(Haddad  et  al.  1994b)  to  develop  complex  structured  singular  value  synthesis  (CSSV) 
formalisms  that  a  priori  fix  the  structure  of  both  the  D-scales  and  the  controller.  Hence, 
an  extension  of  the  algorithms  here  will  enable  fixed-structure  CSSV  controller  synthesis 
that  blends  and  H°°  performance  objectives. 

Practical  applications  often  lead  to  large  dense  systems  of  nonlinear  equations  which 
are  time-consuming  to  solve  on  a  serial  computer.  For  these  systems,  parallel  processing 
may  be  the  only  feasible  means  to  achieving  solution  algorithms  with  acceptable  speed. 
One  economical  way  of  achieving  parallelism  is  to  utilize  the  aggregate  power  of  a  network 
of  heterogeneous  serial  computers.  In  industrial  environments  where  interactive  design  is 
often  the  practice,  the  parallel  code  can  be  easily  incorporated  into  interactive  software 
such  as  MATLAB  or  Mathematica  with  proper  setup  of  the  network  computers.  To  the 
engineering  users  the  design  environment  is  identical.  However  the  computations  are  faster. 

The  most  expensive  part  of  the  homotopy  algorithm  of  (Ge  et  al.  1994)  is 

the  computation  of  the  Jacobian  matrix,  which  can  be  parallelized  easily  to  run  across 
an  Ethernet  network  with  little  modification  of  the  original  sequential  code,  and  which 
has  relatively  large  task  granularity.  There  is  a  trade-off  between  the  programming  effort 
and  the  speedup  of  the  parallel  program.  To  obtain  a  better  speedup,  other  parts  of  the 
homotopy  algorithm,  such  as  finding  the  solution  to  the  Riccati  equations  and  the  QR 
factorization  to  compute  the  kernel  of  the  Jacobian  matrix,  need  to  be  parallelized  as  well. 

In  this  study  the  homotopy  algorithm  for  controller  synthesis  is  parallelized  to 

run  on  a  network  of  workstations  using  PVM  (Parallel  Virtual  Machine)  and  on  an  Intel 
Paragon  parallel  computer,  under  the  philosophy  that  as  few  changes  as  possible  are  to  be 
made  to  the  sequential  code  while  achieving  an  acceptable  speedup.  The  parallelized  com¬ 
putation  is  that  of  the  Jacobian  matrix,  which  is  carried  out  in  the  master-slave  paradigm 
by  functional  parallelism,  that  is,  each  machine  computes  a  different  column  of  the  Jaco¬ 
bian  matrix  with  its  own  data.  Unless  the  Riccati  equation  solver  is  parallelized,  there  is  a 
large  amount  of  data  needed  for  each  slave  process  at  each  step  of  the  homotopy  algorithm. 
To  avoid  sending  too  many  large  messages  across  the  network  or  among  different  nodes  on 
the  Intel  Paragon,  all  slave  processes  repeat  part  of  the  computation  done  by  the  master 
process,  which  therefore  decreases  the  speedup  of  the  parallel  computation. 

The  speedups  of  the  parallel  code  are  compared  as  the  number  of  workstations  increases 
or  as  the  number  of  nodes  increases  on  an  Intel  Paragon  and  as  the  size  of  the  problem  varies. 
A  reasonable  speedup  can  be  achieved  using  an  existing  network  of  workstations  compared 
to  that  of  using  an  expensive  parallel  machine,  the  Intel  Paragon.  It  is  demonstrated 
that  for  a  large  problem,  the  approach  of  using  a  network  of  workstations  to  parallelism  is 
feasible  and  practical,  and  provides  an  efficient  and  economical  computational  method  to 
parallelize  a  homotopy  based  algorithm  for  controller  synthesis  in  a  workstation- 

based  interactive  design  environment. 
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-  (1) 


(2) 


2  The  controller  synthesis  problem 

The  LQG  controller  synthesis  problem  with  an  H°°  performance  bound  can  be  stated  as 
follows:  given  the  ?i-th  order  stabilizable  and  detectable  plant 

x{t)  =  Ax(t)  + Bu{t)  +  Diw{t), 
y{t)  =  Cx{t)  +D2w{t), 

where  A  €  B  G  C  €  R'^”,  £>i  €  €  R'^p,  DiD'^  =  0,  and  w{t)  is 

p-dimensional  white  noise,  find  a  Uc-th  order  dynamic  compensator 

Xc{t)  =  Ac  Xc{t)  +  Be  y{t), 

CcXc{i‘)i 

where  Ac  €  Be  €  Cc  €  and  iic  <  n,  which  satisfies  the  following 

criteria: 

(i)  the  closed-loop  system  (l)-(2)  is  asymptotically  stable,  i.e.,  A  —  ^A^^  ^ 

asymptotically  stable; 

(ii)  the  Qoo  X  P  closed-loop  transfer  function 

H{s)  =  Ec^{sh-A)-^b, 

from  w{t)  to 

z{t)  -  EicoX{t)  -I-  E2ooU{t),  (3) 

where  =  ( ^loo  E2ooCc )  (jBioo  €  R«~^”,  jEjoc  €  £^200  =  0),  n  =  n  +  Uc, 

,  satisfies  the  constraint 

l|ffWlloo<7  (4) 

where  7  >  0  is  a  given  constant;  and 

(iii)  the  performance  functional 

J(Ac,Bc,Cc)  =  lim  £  \x^(t)Rix{t)  -1-  u'^{t)R2u{t)]  (5) 

t—*-00  *’ 

is  minimized,  where  £  is  the  expected  value,  Ri  =  EfEi  €  R"’^”  and  R2  =  E2E2  €  R’"’^’” 
(El  G  R’^",  E2  €  R'*'^"*,  E[E2  =  0  )  are,  respectively,  symmetric  positive  semidefinite 
and  symmetric  positive  definite  weighting  matrices. 

The  closed-loop  system  (l)-(3)  can  be  written  as  the  augmented  system 

x{t)  =  Ax{t)  +  bw{t), 

z(t)  =  EcoX{t) 


(6) 


where  x 


-  (;)■ 


Using  this  notation  and  under  the  condition  that  A  is  asymptotically  stable,  for  a  given 
compensator  the  performance  (5)  is  given  by 

J{Ac,Bc,Cc)  =  ir{QR),  f7) 

where  R=  C^R  C  )  ^  satisfies  the  Lyapunov  equation 
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with  symmetric  positive  semidefinite  Vi  =  DiDf,  symmetric  positive  definite  V2  =  , 

and 

^=(0  bAbj)- 

Lemma  1  (Bernstein  and  Haddad  1989),  (Haddad  and  Bernstein  1990). _ Let 
(Ac-BcjCc)  be  given  and  assume  there  exists  Q  €  satisfying 

Q  is  symmetric  and  nonnegative  definite  (9) 

and 

AQ  +  +  7”^  QRoo  Q  +  V  =  0,  (10) 

where  Roo  =  C  )  ’  ^  and  R^oo  =  -^^00^200  are  symmetric 

positive  semidefinite  matrices.  Then 

(A,jD)  is  stabilizable  (11) 

if  and  only  if 

A  is  asymptotically  stable. 

In  this  case 

||i?(5)||oo  <  7,  (12) 

Q  ^  Q  (Q  —  Qis  nonnegative  definite), 

and 

tr  QR  =  J{Ac,  Be.  Ce)  <  J{Ae,  Cc)  =  tr  QR. 


Hence  the  satisfaction  of  (9)  and  (10)  along  with  the  generic  condition  (11)  leads  to:  1) 
closed-loop  stability;  2)  prespecified  H°°  attenuation;  and  3)  an  upper  bound  for  the 
performance  criterion. 

#  The  auxiliary  minimization  problem  is  to  determine  (Ac,  Be,  Ce)  that  minimizes 

J{Ae,  Be,  Ce),  and  thus  provides  a  bound  for  the  actual  criterion  J(Ac,  Be,  Ce). 

For  technical  reasons  {Ac,Be,Ce,  Q)  is  restricted  to  the  open  set 

S  =  ((Ac,Hc,  Gc,  Q) :  A  and  A  -I-  7  ^QR  are  asymptotically  stable, 

A  6  is  symmetric  positive  definite,  and  {Ae,Be,Cc)  is  controllable  and  observable  }. 


3  The  homotopy  algorithm 


Treating  all  the  components  of  Ac,  Be,  Ce  as  unknowns  is  redundant,  since  only  rieim  +  l) 
parameters  suffice  to  describe  the  full  or  reduced  order  controller.  There  are  numerous 
ways  to  parametrize  the  model  (Ac,  Be,  Ce)  —  the  parametrization  introduced  by  (Ly  et 
al.  1985)  is  one  way  that  uses  the  minimal  number  of  parameters.  In  this  basis,  Ac  is  a 
2x2  block-diagonal  matrix  (2x2  blocks  with  an  additional  1x1  block  jf  Ue  is  od4)  with 


2x2  blocks  in  the  form 


Be  is  a  full  matrix,  and  Gc  has  its  first  row  fixed  as 


(10  1  0  •  •  • ).  Observe  that  the  Ly  structure  has  nc(m  + 1)  unknowns  —  rie  from.A-c) 

Uel  from  Be,  and  nc(m  -  1)  from  Ce-  Denote  the  unknowns  by 


/  iAc)j  \ 

0  =  Vec  (Be)  , 
V  Vec  (Gc)r.  / 
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where  {Ac)j  is  a  vector  consisting  of  those  elements  in  Ac  indexed  by  the  set  2  = 
{(2, 1),  (2, 2), . . . ,  (ric,  Uc)},  Vec  of  a  matrix  is  the  concatenation  of  its  columns,  and  (C'c)r- 
is  the  matrix  obtained  from  rows  T  =  {2, . . .  ,m}  of  Cc- 

Choose  a  problem  for  which  a  solution  is  known,  defined  by  the  matrices  j4o,  Bq, 
Co,  i^i.o,  ^2,0,  -Rioo.o,  ■R2oo,0,  14,0,  14.0,  7o;  sxactly  how  this  initial  problem  is  chosen  is 
described  in  (Ge  et  al.  1994).  Let  Af,  Bf,  ...,  7/,  denote  A,  B,  7,  in  the  above  and 
define  >1(A),  . . .,  7(A),  as 

A{X)  =  Ao  +  X{Af  -  Ao),  7(4)  =  70  +  4(7/ -  7o), 

and  for  brevity  denote  them  by  A,  . . .,  7,  respectively  in  the  following.  Let 
HAc{e,\)  =  Vl2Ql2+'P2Q2. 

HbAO,  4)  =  V2BCV2  +  {TI2Q1  +  V2Qj2)C^, 

HcS^i  4)  =  R2CCQ2  +  B^{'PiQi2  +  ^1262) 

+  ^7"'i?2ooC'4(Qf2Pi  -f-  Q2Vl2)Ql2  +  (Qr2^12  +  Q2^2)S2], 
where  Q  and  V  satisfy 

AQ+QA^ +  -f-^QRooQ  +  V  =  0,  (13) 

(i  +  i-'^QRoofv  +  V{A  -b  j-^QR^)  +R  =  0,  (14) 


Q  = 


Qi2\ 

Q2J’ 


V12 

V2 


Roo^^  q  C^i?2ooCcj’  Rloo  =  EfooRioo,  R200  =  EJ^E2oo- 

The  nonlinear  equations  (corresponding  to  A  =  1)  that  result  from  the  nonconvex 
constrained  minimization  procedure  (see  (Ge  et  al.  1994)  for  details)  are 

/  \ 
p(0,A)=  Vec  [ifB,(^,A)]  =0. 

VVec  [ffc.(^,4)]^.)/ 

To  guarantee  a  full  rank  Jacobian  matrix  along  the  whole  homotopy  zero  curve  except 
possibly  at  the  solution  corresponding  to  A  =  1,  define  the  homotopy  map  to  be  p{6,  A)  = 
Xp{0,  A)  +  (1  —  A)(0  —  do).  The  Jacobian  matrix  of  p  is  given  by 

Dp{e,  X)  =  {XDep{e,  A)  +  (1  -  A)/,  p{e,  X)  +  XDxpie,  x)-{e-  do)) . 

Dop{d,  A)  and  D\p{d,  A)  are  derived  from  the  following.  Define 

Hac  =  'Pl2^^^Qi2  +  Vl2Qn  +  +  V2Q^^\ 

=  v\^^BcV2  +  +  v\'''Ql2  +  V2Qli'^)C'^, 

=  R2CcQ^i^  + 

+  ^i22ooCe[(e[y^Pl  +  +  QfRn  +  Q2V'[2^'^)Qi2  +  {Qn^^'^ru  +  Ql2'Pl2 

+  QfV2  +  Q2V^^)Q2  +  (Qfa^l  +  02pr2)ei2  +  {QnRn  +  S2:P2)Q^2'^], 
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where  the  superscript  (j)  means  djdOj.  Using  the  above  definitions,  we  have  for  6j 
(ylc)fep  where  {k,l)  €  I, 


d{A,)u 


for  Oj  =  {Bc)^ 


0{Ac)kl 

0{^c)kl 


=  Hb^V^^K  +  'P2E^'‘''^V2. 

d[Bc)kl 


and  for  $j  =  {Cc)kh  where  fc  >  1, 


diC^h 


where 


v(C/c)fc^ 

+  R2E^*^''-^Q2  +  ^R2ooE<^'^'‘^Y, 
o{Cc)kl  ^ 


Y  =  {Ql^Vl  +  Q2PI2)  2i2  +  (Qr2^12  +  Q2^2)Q2 


and  is  a  matrix  of  the  appropriate  dimension  whose  only  nonzero  element  is  tki  —  1. 
and  Oj'A  can  be  obtained  by  solving  the  Lyapunov  equation 

0  =  (i  +  ')-^QRj)  +  Q^^\A  +  'y-^QRoof  + 

+  Q  +  +  'r~^QRg^Q, 

0  =  (i  +  'r-^QRoofv^'>  +  v^'>  {A  +  j-^QRoo)  +  R^^'> 

+  (A<J>  +  'r-^Q^^'>Roo  +  'y-^QR^^fv  +  V{A^^'^  + 

Similarly  for  A,  using  a  dot  to  denote  9/5A, 

^  =  Ha.(A2), 

^  =  Hb^  {V,  Q)  +  V2B,V2  +  {Vj^  Qi  +  V2  QJ2)C^^ 

=  iicciPi  2)  +  R2CCQ2  +  {R1Q12  +  R12Q2)  +  2^  ^R2oo('cY—'y  ^jR2ooCcY, 
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where  V  and  Q  are  obtained  by  solving 

.  .  .  ^  - 

0  =  (a -f- 7~^Q^oo) Q  +  Q(-d -b 7  ^QRoo)  +V  +  AQ+QA  -I-7  "^QRooQ 

-  2-y~^'y QRooQ, 

0  =  (A  -t-  'Y'^QRoofv  +  V{A  +  'y-^QRoo)  +  ^  +  (A  +  -f~^QRoo  +  ' 

-  2y-^QRoofv  +  V{A  +  'y-^QRoc  +  -  2j-^QRoo)- 

The  homotopy  zero  curve  tracking  algorithm  (which  is  a  standard  globally  convergent 
probability-one  homotopy  algorithm  (Watson  et  at  1987))  is 

1)  Set  A  :=  0,  0  :=  ^o- 

2)  Calculate  R,  Roo,  V,  and  compute  Q  and  V  according  to  (13)  and  (14). 

3)  Evaluate  the  homotopy  map  p(0,  A)  and  the  Jacobian  of  the  homotopy  map  Dp{6,X). 

4)  Predict  the  next  point  on  the  homotopy  zero  curve  using,  e.g.,  a 

Hermite  cubic  interpolant. 

5)  For  fc  :=  0, 1, 2,  •  •  •  until  convergence  do 

^(fe+i)  ^  ^(fe)  _  [Dp{Z^'‘'>)]^ p{Z^’^'>), 

where  [jDp(Z)]^  is  the  Moore-Penrose  inverse  of  Dp{Z).  Let  (^i,Ai)  =  Jim 

6)  If  Ai  <  1,  then  set  6  :=  61,  X  :=  Ai,  and  go  to  step  2). 

7)  If  Ai  >  1,  compute  the  solution  0  at  A  =  1. 

4  The  parallel  algorithm 

In  the  above  algorithm,  Step  2)  involves  solving  one  Riccati  equation  and  one  Lyapunov 
equation.  The  Riccati  equation  is  solved  using  Laub’s  Schur  method  (Laub  1979).  The 
algorithm  of  Bartels  and  Stewart  (Bartels  and  Stewart  1972)  is  applied  to  solve  the  Lya¬ 
punov  equation.  Although  both  algorithms  are  0((n  +  Ucf),  the  Riccati  equation,  being 
more  complicated,  takes  much  more  CPU  time  to  solve.  Once  Q  and  V  are  obtained,  the 
homotopy  map  p  is  formed  by  matrix  multiplication  operations. 

The  major  part  of  the  computation  in  Step  3)  is  that  of  the  Jacobian  matrix.  The 
number  of  variables  including  A  in  this  formulation  is  ndw,  + 1)  -f  1.  Each  column  of  the 
Jaoobian  matrix  corresponds  to  the  derivative  of  the  homotopy  map  with  respect  to  one 
variable  and  requires  the  solution  of  two  Lyapunov  equations  (Ge  et  at  1994).  Therefore 
the  time  complexity  of  the  Jacobian  matrix  computation  is  0{nc{rn  4-  l){n  -1-  tIc)^).  The 
Bartels  and  Stewart  algorithm  finds  the  real  Schur  form  of  A  or  A^  depending  on  the 
Lyapunov  equation.  At  each  step  along  the  homotopy  path  unnecessary  factorization  can 
be  avoided  if  the  previous  factorization  results  from  the  computation  of  p  and  Dxp  are 
used. 

The  primary  goal  of  this  study  is  to  make  use  of  the  existing  code  and  to  achieve  rea¬ 
sonable  parallel  eflBlciency  economically.  The  only  part  of  the  algorithm  that  is  parallelized 
is  the  Jacobian  matrix  computation  in  Step  3).  To  utilize  existing  computer  resources  such 
as  a  network  of  workstations,  the  software  package  PVM  (Parallel  Virtual  Machine)  is  used 
to  provide  the  distributed  computing  capabilities. 

The  parallel  algorithm  follows  the  master-slave  paradigm.  The  master  sends  the  index 
of  the  column  of  the  Jacobian  matrix  to  be  computed  to  a  slave.  The  slave  computes  the 
corresponding  column  of  the  Jacobian,  sends  the  column  back  to  the  master,  and  waits  for 
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Number  of  worstations  Number  of  worstations 

Fig.  1.  Speedup  with  master  and  slaves  Fig.  2.  Speedup  when  one  slave  is  on 
on  different  machines.  the  master  machine. 

the  next  index  from  the  master  to  arrive.  After  receiving  a  column  of  the  Jacobian,  the 
master  sends  another  index  to  the  idle  slave.  In  the  implementation  for  the  Intel  Paragon, 
asynchronous  send  which  sends  a  message  without  waiting  for  completion  is  used  whenever 
possible  to  speed  up  the  communication. 

When  the  algorithm  is  implemented  on  a  network  of  workstations,  the  modification  to 
the  original  sequential  source  code  consists  of  three  parts:  the  first  one  is  to  spawn  slave 
processes  and  set  up  the  communication  links  between  the  master  and  the  slaves;  the  second 
is  to  extract  a  slave  program  from  the  original  code  and  at  the  same  time  simplify  the  master 
program;  the  last  is  to  add  a  mechanism  to  guarantee  correct  communication  between 
master  and  slaves.  The  first  part  consists  of  standard  PVM  operations,  while  the  second  is 
more  problem  oriented.  To  decrease  communication,  each  slave  process  repeats  part  of  the 
computation  of  p  and  Dxp  so  that  Q  and  V  are  not  sent  through  the  network.  There  is  no 
loss  of  eflBciency  since  the  master  is  also  computing  the  same  quantities.  The  slave  program 
consists  of  mainly  the  original  subroutines  with  additional  code  for  communication. 

For  the  implementation  on  the  Intel  Paragon,  the  modification  of  the  original  code 
is  even  simpler.  There  is  no  need  for  a  separate  slave  program  if  control  statements  use 
node  identification  properly.  The  parent  process  run  on  an  Intel  Paragon  always  gets  node 
number  0,  while  other  nodes  are  numbered  1  and  higher.  The  statement  if  node-Dumber 
==  0  precedes  the  code  that  is  to  be  executed  by  the  master,  and  an  else  following  the 
previous  master  code  will  precede  the  code  to  be  executed  by  the  slave.  The  remaining 
modification  to  the  original  code  is  similar  to  the  implementation  using  PVM.  Asynchronous 
send  is  used  whenever  possible.  A  wait  is  used  later  when  the  data  is  needed,  to  ensure 
correct  communication  between  the  master  and  the  slaves. 

5  Results 

The  distributed  code  using  PVM  was  run  on  a  network  of  seven  SGI  Indigo^  worksta¬ 
tions.  The  data  came  firom  a  control  problem  for  suppressing  vibrations  in  a  string  under 
transverse  loading  firom  a  time  varying  disturbance  force  (Erwin  1993).  For  dimensions 
71  =  12,  20,  28,  36  reduced  order  controllers  of  dimensions  Uc  =  10,  18,  26,  34  are  sought 
respectively. 

The  speedups  versus  the  number  of  workstations  are  shown  in  Figs.  1  and  2  {n  =-36, 
28,  20,  12,  top  to  bottom).  Fig.  1  shows  the  speedup  versus  the  number  of  workstations 
when  the  master  process  and  the  slave  processes  are  run  on  different  machines,  while  Fig.  2 
corresponds  to  the  situation  where  the  master  process  and  a  slave  process  with  lower  priority 
are  run  on  one  machine  and  the  rest  of  the  slaves  are  on  other  machines.  For  fair  comparison 
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Fig.  3.  Comparison  of  speedups. 

all  the  speedups  are  computed  relative  to  the  results  of  the  best  optimized  sequential  code. 
In  Fig.  1  two  workstations  correspond  to  the  master  process  on  one  machine  and  the  only 
slave  on  the  other. 

As  shown  in  the  figures,  the  speedup  increases  as  the  dimension  of  the  problem  increases, 
or  as  the  number  of  workstations  increases  for  a  sufficiently  large  problem.  The  speedups 
from  three  scenarios  (solid  line — master  and  slaves  on  different  machines;  dash-dot  line — 
master  and  a  low  priority  slave  on  one  machine  and  the  rest  of  the  slaves  on  others;  dashed 
line — master  and  a  slave  with  the  same  priority  as  the  master  on  one  machine)  for  n  =  20 
are  plotted  against  the  number  of  workstations  in  Fig.  3.  If  the  number  of  workstations 
is  <  4  it  is  better  to  use  the  second  scenario.  When  the  number  of  workstations  is  >  4, 
the  speedup  is  higher  if  all  the  processes  including  the  master  and  the  slaves  are  run  on 
different  machines.  Similar  results  obtain  for  large  n. 

The  same  algorithm  is  implemented  using  the  system  function  calls  of  an  Intel  Paragon 
and  run  on  one  with  28  processors  at  Virginia  Polytechnic  Institute  and  State  University. 
Fig.  4  shows  the  results  obtained  from  the  run  for  n  =  12,  20,  28.  The  number  of  nodes 
varies  up  to  25.  The  highest  curve  corresponds  to  the  speedup  when  n  =  28  and  the  lowest 
corresponds  to  that  when  n  =  12.  As  n  increases,  the  advantage  of  parallel  processing  also 
increases.  The  highest  speedup  achieved  for  n  =  28  using  seven  SGI  Indigo^  workstations 
is  about  3.3  while  the  highest  speedup  using  25  nodes  on  an  Intel  Paragon  is  about  5.1. 
Comparing  speedups  is  meaningful  since  the  performance  of  a  single  SGI  Indigo^  processor 
is  roughly  comparable  to  that  of  a  single  i860XP  Paragon  node  (actually,  depending  on  the 
task,  the  lOOMHz  R4000  Indigo^  is  faster  by  a  factor  of  2).  However,  the  cost  of  the  Intel 
paragon  is  a  factor  of  three  times  the  cost  of  the  SGI  UNIX  workstation  network.  Much 
higher  speedups  are  potentially  possible  with  the  Paragon,  but  not  without  considerable 
programming  effort  for  this  controller  design  problem. 

The  above  methodology  can  be  easily  generalized  to  industrial  design  environments 
where  software  packages  like  MATLAB  or  Mathematica  are  often  used.  The  sequential 
program  for  mixed-norm  LQG  controller  synthesis  has  been  developed  into  a 

MATLAB  package.  It  is  easy  to  include  this  distributed  implementation  into  the  MATLAB 
package.  Installation  requires  two  steps:  the  first  one  is  to  install  PVM  on  the  network 
of  workstations,  and  the  second  is  to  create  a  file  in  which  all  the  worker  machines  on 
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Fig.  4.  Results  for  Intel  Paragon  XPE-28. 

the  network  are  listed  (Geist  et  al.  1993).  The  execution  of  the  distributed  program  from 
within  an  interactive  design  environment,  e.g.,  MATLAB,  can  be  done  by  using  a  MATLAB 
function  defined  in  a  MATLAB  .m  file,  in  which  UNIX  shell  commands  will  start  the  PVM 
daemons  if  they  have  not  already  been  started,  and  will  then  execute  the  distributed  code. 

In  summary,  parallelizing  the  Jacobian  matrix  computation  in  a  homotopy  algorithm 
reduces  the  execution  tiime  and  is  economical,  especially  for  large  problems.  Acceptable 
speedups  are  obtained  for  a  PVM  implementation  on  a  network  of  workstations.  The 
approach  can  be  applied  to  real  industrial  design  environments  to  reduce  controller  design 
time  and  efl:'ectively  utilize  existing  workstation  networks.  Compared  to  the  cost  of  using 
real  parallel  computers  and  the  cost  of  developing  highly  efficient  parallel  code,  the  approach 
of  utilizing  a  network  of  workstations  with  only  moderately  efficient  code  is  much  more  cost- 
effective  for  controller  synthesis. 
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An  Object-oriented  Approach  to  Semidefinite  Programming 

Yuzhen  Ge+  Layne  T.  Watson*  Emmanuel  G.  Collins,  Jr.* 


Abstract 

An  object-oriented  design  and  implementation  of  a  primal-dual  algorithm  for  solv¬ 
ing  the  semidefinite  programming  problem  is  presented.  The  advantages  of  applying 
the  object-oriented  methodology  to  numerical  computations,  in  particular  to  an  inte¬ 
rior  point  algorithm  for  semidefinite  programming,  or  for  solving  other  types  of  linear 
matrix  inequalities  are  discussed.  One  object-oriented  design  of  the  primal-dual  al¬ 
gorithm  and  its  implementation  using  C-f-H-  is  presented.  The  performance  of  the 
C-h+  implementation  is  compared  with  that  of  a  procedural  C  implementation,  and 
while  the  performance  of  the  0+- f  implementation  is  comparable  to  that  of  the  C 
implementation,  the  resulting  code  is  easier  to  read,  modify,  and  maintain. 

Key  words:  linear  matrix  inequality,  object  oriented,  primal  dual  algorithm,  scientific  com¬ 
putation,  semidefinite  programming 

CR  categories:  D.2.2,  G.1.6 
1  Introduction 

Object-oriented  design  and  programming  has  been  a  major  theme  in  software  engi¬ 
neering  in  recent  years.  Traditional  design  or  classical  design,  which  has  been  the  main 
software  design  paradigm  until  about  mid  80’s,  concentrates  on  the  actions  that  a  system 
has  to  take  and  decomposes  the  system  into  separate  units  or  modules  according  to  their 
functionalities.  In  object-oriented  design  a  system  to  be  modeled  is  viewed  as  a  collec¬ 
tion  of  objects,  each  of  which  has  its  own  attributes  and  the  operations  performed  on  an 
object  or  functions  acting  on  an  object  are  also  defined  in  one  syntactic  unit.  Objects 
communicate  by  passing  messages  or  by  calling  functions  fi:om  other  objects  which  provide 
services.  Object-oriented  design  is  developing  an  object-oriented  model  of  a  system  and  can 
be  realized  (implemented)  by  object-oriented  programming  using  languages  such  as  C-h-b, 
FORTRAN  90,  or  Smalltalk. 

The  advantages  of  object-oriented  design  and  programming  have  been  described  widely 
elsewhere  (!].  A  short  summary  will  be  provided  here.  First,  an  object  is  an  independent 
entity  that  is  encapsulated  in  one  syntactic  unit.  The  definition  of  an  object  consists  of  the 
definition  of  the  attributes  of  the  object  along  with  operations  that  can  be  performed  on 
the  object  and  the  services  or  function  calls  provided  by  the  object.  Encapsulation  hides 
the  implementation  details  of  an  object  and  makes  the  program  easier  to  read  and  modify. 
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Any  subsequent  change  to  the  program  can  be  locahzed,  making  the  resulting  program 
more  easily  maintained. 

The  second  advantage  is  information  hiding.  Definitions  of  an  object  which  need  not 
be  known  to  other  objects  are  inaccessible  to  other  objects,  preventing  them  frona  being 
changed  accidentally.  In  other  words,  information  hiding  makes  implementation  details  of 
an  object  inaccessible  to  other  objects.  However,  the  designer  has  the  freedom  to  decide 
what  to  hide  and  what  not  to  hide. 

The  third  advantage  is  code  reuse.  Inheritance  enables  the  definition  of  a  new  object, 
which  can  be  viewed  as  a  subclass  of  an  existing  object,  without  having  to  repeat  some 
of  the  details.  The  new  object  can  inherit  attributes  or  operations  from  its  ancestor. 
Inheritance  is  one  way  to  support  reuse  of  existing  objects.  There  are  different  kinds  of 
reuse  in  object-oriented  programming;  inheritance  is  only  one  of  them. 

One  of  the  most  popular  object-oriented  programming  languages  is  C-f-+  [11],  which  is 
used  to  implement  the  algorithm  of  this  paper.  Some  of  the  reasons  why  C-I-+  is  so  widely 
used  are  upward  compatibility  with  C,  design  emphasis  on  eSiciency  and  performance,  and 
the  availability  of  many  useful  libraries  and  tools.  For  instance,  the  Gnu  C-H-f  compiler 
and  other  tools  are  available  on  a  wide  range  of  platforms  and  provide  good  performance, 
programming  environments,  and  reasonable  compliance  with  ANSI  standards. 

There  are  many  available  libraries  such  as  IML-t— 1-[6],  SparseLib-l-+  [5]  [9],  STL[10]  [8], 
and  others  which  emphasize  numerical  computation.  One  notable  package  is  LAPACK-b-h, 
developed  by  Dongarra  et  al.  [4],  which  is  a  C-t-f  interface  to  LAPACK  and  BLAS.  Ref. 
[4]  has  shown  that  performance  of  programs  using  the  package  is  comparable  to  calling 
LAPACK  and  BLAS  directly,  and  can  at  the  same  time  reap  the  benefits  of  object-oriented 
programming. 

This  paper  contains  the  result  of  object-oriented  design  and  implementation  of  an  al¬ 
gorithm  for  semidefinite  programming.  Semidefinite  programming  refers  to  minimizing  a 
linear  function  subject  to  a  linear  matrix  inequality  [12],  That  is, 

minimize  (Fx 
subject  to  F{x)  >  0, 

where 

m 

F{x)  =Fo  +  '^XiFi, 

»=i 

c  e  R”*,  and  Fq, , Fm  €  R”^”  are  symmetric  matrices.  F(x)  >  0  means  that  F(x)  is 
positive  semidefinite. 

Many  problems  in  controls  engineering  can  be  cast  in  terms  of  a  semidefinite  program¬ 
ming  problem  [12].  Since  a  semidefinite  programming  problem  is  a  convex  oi)timization 
problem,  which  can  be  solved  by  interior  point  methods  [7],  it  has  attracted  the  attention 
of  many  researchers  in  interior  point  methods.  There  is  a  C  implementation  of  a  primal-dual 
algorithm  for  solving  the  semidefinite  programming  problem  [13].  A  C-1-+  implementation 
of  that  primal-dual  algorithm  for  the  semidefinite  programming  problem  is  developed  here 
to  explore  the  possible  benefits  of  object-oriented  design.  Because  of  the  similarity  of  the 
primal-dual  algorithm  with  other  interior  point  algorithms  for  solving  the  semidefinite  pro¬ 
gramming  problem,  the  design  and  implementation  methodology  developed  here  can  be 
easily  modified  and  applied  to  other  interior  point  algorithms. 
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The  performance  of  a  C-f-+  implementation  of  the  primal-dual  algorithm  for  semidefinite 
programming  is  compared  with  the  existing  C  implementation  of  the  same  algorithm  firom 
(13].  While  the  CPU  times  of  the  two  implementations  are  comparable  to  each  other,  the 
C++  version  oflfers  the  advantages  mentioned  earlier  in  this  section.  Segments  of  the  code 
will  be  used  to  illustrate  object-oriented  features  of  the  implementation. 

Section  2  briefly  sketches  the  primal-dual  algorithm  for  semidefinite  programming  that 
will  be  used  to  illustrate  the  object-oriented  design  methodology.  In  Section  3,  the  details  of 
an  object-oriented  design  of  the  primal-dual  algorithm  will  be  given.  The  implementation 
using  C++  will  be  described  in  Section  4.  Comparison  and  discussion  of  the  C  and  C++ 
results  will  be  given  in  Section  5. 


2  Primal-dual  algorithm  for  semidefinite  programming 

The  primal-dual  algorithm  for  solving  the  semidefinite  programming  program  given 
in  detail  in  [12]  will  be  described  briefly  here.  The  dual  problem  associated  with  the 
semidefinite  program  (1)  is  ^ 

maximize  —\,xFoZ 

subject  to  tr  FiZ  =  Cj,  z  =  1, . . . , m  (2) 

Z  =  Z  >  0. 


The  primal-dual  method  can  be  interpreted  as  solving  the  primal-dual  optimization 
problem 

minimize  c^x  +  tTFoZ 
xeR'",Z€R’*’‘" 

subject  to  tv  FiZ  =  Ci,  z  =  1, . . . ,  m 
F{x)  >0,Z>0,Z  =  Z^, 


where  the  objective  function  c^x  +  tr  FqZ  =  z;  is  called  the  duality  gap,  which  has  the 
known  optimum  value  of  zero  for  a  convex  problem.  The  advantage  of  using  the  primal- 
dual  formulation  is  that  at  each  step  information  from  the  dual  problem  can  be  used  to 
obtain  a  good  update  for  the  primal  variables. 

One  of  the  methods  to  solve  the  primal-dual  optimization  problem  is  the  potential 
reduction  method.  Define  a  potential  function 

(f>(x,  Z)  =  (n  +  vi/n)  log(  tr  F{x)Z)  —  logdet  F(x)  —  log  det  Z  —  nlogn, 


where  f  >  1  is  a  weighting  parameter  in  the  potential.  The  duality  gap  is 


rj  <  exp 


which  is  small  if  the  potential  function  (f>  is  small  and  approaches  0  if  is  approaching  — oo. 

The  whole  algorithmic  process  can  be  described  as  follows.  Starting  from  a  strictly 
feasible  xo  and  Zq,  find  Xk  and  Zfc  so  that  the  potential  <t>  is  reduced  at  each  step  by  at 
least  a  fixed  amount  6  in  every  step 


until  the  duality  gap  rj  is  smaller  than  some  specified  e  >  0.  The  first  x*  and  Zfc  which 
make  z;  <  €  is  the  approximate  numerical  solution  of  the  primal-dual  problem. 


SpAlgo 
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1 1  SpPrime  prime 
i  (  SplHiel  «*dual 
number  p«  q 
cTx<  3 •  traceZF0( ) $ 
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SearchRectC • . « ) 
updateMar ( . . < ) 
PlaneSearch< . • . ) 


Gapp 

.)  I 


I S  prine.variable  ^ 

' {  matrix  c  i 
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I  xnormf  xiiom«  norm* 
traf»ZFB<3«  gapCx)  | 
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Searchltect(p«  q«  « . .  3  \ 
Plansearch(p«  q««..3  i 
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iS  dual_variable  Z  ^ 
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f  a  sum_Fi_x  * 

traceZFB(3a  gap<x3  ' 

ReducePotential ( . . . 3  \ 

SearchRecrtCpf  qa>--3  ^ 

PlansearchCp#  qa...3  | 


SpSparse? 


Fig.  1.  Class  diagram. 

The  primal-dual  potential  reduction  method  for  solving  the  semidefinite  program  (1) 
can  be  summarized  as  follows.  Given  strictly  feasible  x  and  Z,  while  duality  gap  = 
c^x  -F  tr  FqZ  >  e  do 

1.  compute  a  direction  6x  for  the  primal  variable  x  and  a  direction  6Z  for  the  dual 

variable  Z;  ,,  , 

2.  find  p,  q  that  minimize  (f){x  +  p6x,  Z  +  q6Z),  where  (p,  q)  are  constrained  to  some 
search  rectangle  in  the  plane; 

3.  update  x  :=x+  p6x  and  Z  :=  Z  q6Z.  ' 

Details  of  each  of  these  steps  is  given  in  [12]. 

3  The  object-oriented  design 

In  C-F- 1-  terminology,  the  word  “class”  is  used  to  mean  a  type  of  object.  For  example,  a 
class  Symm  can  be  defined  for  symmetric  matrices,  while  a  specific  symmetric  matrix  is  an 
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object  of  type  Symm  or  an  instance  of  the  class  Symm.  C-f-+  terminology  will  be  followed 
in  the  rest  of  the  paper. 

Object-oriented  analysis  and  design  is  one  of  the  most  active  research  areas  for  both 
academics  and  industrial  practitioners.  When  applied  in  different  circumstances,  different 
analysis  and  design  techniques  may  be  emphasized.  One  of  the  most  influential  analysis 
and  design  techniques  is  due  to  Boochfl],  whose  conventions  will  be  loosely  followed  in  this 
paper. 

First  we  will  describe  the  classes  used  in  the  design  and  their  relationship.  The  class 
diagram  for  the  problem  is  shown  in  Fig.  1,  where  all  the  classes  are  defined  with  their 
attributes  and  functions.  For  clarity,  only  main  attributes  and  functions  are  shown  and  the 
functions  names  may  be  different  from  that  used  in  real  implementation. 

In  Fig.  1,  each  dotted  cloud-shaped  figure  describes  a  class,  with  the  private  members 
of  the  class  preceded  by  ||  and  the  protected  members  of  the  class  preceded  by  |.  Protected 
properties  are  inherited  by  and  accessible  to  subclasses,  whereas  private  properties  are  not. 
AH  other  properties  are  considered  public,  i.e.,  accessible  by  other  classes.  Protected  and 
private  properties  are  not  accessible  to  other  classes.  For  example,  in  the  class  SpPrime,  the 
private  attributes  are  prime-variable  x,  matrix  c,  and  increment  dx,  the  main  public  func¬ 
tions  are  cTx()  which  computes  c^x,  cTdxQ  which  computes  c^cte,  Dx()  which  computes 
dx,  and  update()  to  update  x. 

The  connection  • - ■  from  SpAlgo  to  SpPrime  denotes  the  physical  containment  of 

SpPrime  in  SpAlgo,  while  the  connection  • — o  from  SpAlgo  to  SpDual  denotes  a  pointer 
reference  to  the  class  SpDual  by  SpAlgo,  which  can  be  illustrated  by  the  corresponding 
part  of  the  definition  of  SpAlgo 


private : 

SpPrime  prime;  //  primal  space 

SpDual  *dual;  //  dual  space,  dynamically  allocated 

This  code  declares  prime  as  an  instance  of  the  class  SpPrime,  and  dual  as  an  instance 
of  the  class  SpDual. 

An  arrow  —*  from  SpSym  to  SpDual  denotes  that  SpSym  is  a  subclass  of  SpDual  and 
inherits  the  public  and  protected  attributes  and  functions  from  SpDual.  An  upside-down 
triangle  with  A  in  SpDual  denotes  that  SpDual  is  an  abstract  class.  Many  functions  in 
SpDual  are  defined  as  virtual  so  that  dynamic  binding  is  used  for  these  functions,  i.e.  the 
decision  on  which  implementation  to  use  is  made  at  run  time,  it  it  cannot  be  determined  at 
compile  time.  Consequently,  if  a  class  other  than  SpSym  is  used  to  implement  SpDual,  then 
the  run  time  system  will  choose  the  right  function  depending  on  the  parameters  passed. 

The  symmetry  between  the  primal  and  dual  spaces  in  the  objective  function  in  (3) 
suggests  that  these  two  spaces  should  be  fundamental  classes  in  the  problem,  and  so  two 
classes,  SpPrime  and  SpDual,  are  defined.  SpPrime  is  relatively  simple.  The  primal  variable 
x  is  the  independent  variable.  The  only  important  actions  on  the  class  are  to  calculate  c^x 
and  update  x.  The  class  SpDual  will  contain  more  functions  and  is  more  complicated.  It  will 
not  only  calculate  tcFaZ,  but  also  compute  6Z,  p,  q,  and  update  Z.  The  dual  variatiTe  Z  is 
a  symmetric  matrix.  In  the  future  we  may  want  to  exploit  any  special  structure  the  matrix 
Z  may  have.  For  example,  we  could  exploit  the  sparsity  of  the  matrix  Z  by  defining  and 
using  another  class  SpSparse.  So  SpDual  is  designed  as  an  “envelope”  class  [3]  to  provide 
a  generic  dual  space  interface  and  to  isolate  it  from  the  “letter”  class  that  implements  the 
dual  space  for  a  particular  representation.  As  of  this  writing,  we  implemented  SpSym,  a 


symmetric  matrix  representation  not  exploiting  any  other  special  structure  that  the  matrix 
may  have. 

SpAlgo  is  the  base  class  for  algorithm  implementation.  It  mediates  the  communication 
between  SpPrime  and  SpDual,  executes  the  computational  tasks,  and  changes  the  “states” 
of  the  SpPrime  and  SpDual  objects. 

Finally  main  is  a  utility  class,  a  “free”  program  not  tied  to  any  particular  class,  that 
stores  the  convergence  criterion  and  is  the  driver  for  the  algorithm.  It  is  denoted  with  a 
dotted  cloud  with  a  shadow  in  Fig.  2  and  Fig.  3. 

The  major  advantage  of  this  way  of  thinking  is  that  the  implementation  of  the  dual 
space  is  completely  independent  of  the  program  structiue  and  state  transition.  This  is 
desirable  because  the  dual  space  implementation  is  where  algorithmic  changes  are  likely 
to  occur,  and  this  separation  of  the  interface  and  implementation  localizes  changes  in  the 
program.  The  envelope  class  SpDual  provides  a  clean  and  effective  interface  for  the  dual 
space;  SpSym  is  one  of  the  possible  subclasses  of  SpDual  that  does  nothing  but  implement 
the  specification  of  SpDual.  The  matrix  and  vector  classes  are  from  LAPACK++,  and  are 
not  shown  in  the  diagram. 

There  are  two  execution  scenarios,  whose  outlines  are  shown  in  Figs.  2  and  3,  corre¬ 
sponding  to  different  stages  of  the  primal-dual  algorithm. 

A  small  square  with  F  inside  on  the  border  of  SpDual  and  SpPrime  on  the  line  from 
SpAlgo  denotes  that  SpPrime  and  SpDual  are  part  of  the  client  object,  SpAlgo.  The  two 
square  boxes  with  L  inside  on  the  border  of  SpAlgo  on  the  line  from  main  denote  that 
SpAlgo  is  a  locally  declared  object  in  main.  An  arrow  with  an  empty  circle  indicates  that 
results  or  services  are  returned  to  the  object  pointed  to  by  the  arrow.  An  arrow  denotes 
the  direction  of  a  message  from  one  object  to  another.  The  number  in  front  of  the  message 
denotes  the  execution  sequence  of  the  message. 


In  the  first  scenaxio  (Fig.  2): 

Step  1.  If  a  search  rectangle  exists,  main  sends  a  message  to  SpAlgo  to  do  a  plane  search,  i.e., 
find  p  and  q  such  that  <f){x+p6x,  Z  +  q6Z)  is  minimized,  where  p  and  q  are  constrained 
in  the  search  rectangle.  Otherwise  main  sends  a  message  to  SpAlgo  to  calculate  the 
duality  gap. 

Step  2.  SpAlgo  sends  a  message  to  SpPrime  to  get  the  primal  variable. 

Step  3.  SpAlgo  sends  the  primal  variable  to  SpDual  along  with  a  message  to  calculate  the 
duality  gap. 

Step  4.  main  checks  the  convergence  criterion.  If  the  criterion  is  met,  execution  stops.  Other¬ 
wise,  go  to  the  second  scenario. 

In  the  second  scenario  (Fig.  3): 

Step  1.  main  sends  a  message  to  SpAlgo  to  do  a  potential  reduction  calculation. 

Step  2.  SpAlgo  sends  a  message  to  SpDual,  which  returns  dx,  the  increment  for  the  primal 
variable. 

Step  3.  SpAlgo  sends  dx  to  SpPrime  to  store  it  in  SpPrime. 

Step  4.  main  sends  a  message  to  SpAlgo  to  calculate  the  search  rectangle. 

Step  5.  SpAlgo  sends  a  message  to  SpPrime  to  calculate  c^dx. 

Step  6.  SpAlgo  sends  (Fdx  along  with  message  SearchRect{...)  to  SpDual  to  calculate  and 
return  several  intermediate  variables  that  mark  the  search  rectangle. 

Step  7.  main  checks  the  convergence  criterion  using  the  corners  of  the  search  rectangle.  If  the 
convergence  criterion  is  met,  execution  stops.  Otherwise,  go  to  the  first  scenario. 
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Fig.  2.  Execution  Scenario  1. 


Fig.  3.  Execution  Scenario  2. 

4  The  C-|-+  implementation 

The  program  is  built  upon  the  LAPACK++  vl.O  paxikage,  especially  the  LaVectorDou- 
ble,  LaGenMat,  LaSymmMat  classes,  and  BLAS-I-+  is  used  extensively.  Several  special 
purpose  routines  are  added  to  the  LAPACK++  package  to  accommodate  the  primal-dual 


for  (k=0,  pos4=pos;  k<blck_szs [i] ;  pos4+=blck_szs[i]-k,  k++) 

{ 

seal  =  8igx[k]  ; 

rh.s[po84]  *  (1,0/scal  +  rho*8cal)/sqrt2; 

} 

for  (j«0,  pos=0;  j<m;  j++) 
for  (1=0,  pos2=0;  KL; 

po8  +=  blck_8zsCi]*(blck_szs[i]+l)/2, 
pos2  +=  blck^szs  [i]  ♦blck.szs  [i] ,  i++) 

/♦  compute  V*  ♦  Fj(i)  ♦  V,  store  In  Fsc+pos,  V  is  scaled.*/ 
ciigmcb(2,  blck.szsCi],  F+sz+pos,  R+pos2,  Fsc+pos,  tenq)); 

/*  correct  diagonal  elements  */ 

for  (k=0,  pos4=pos;  k<blck.szs[i] j  pos4  +=  blck_szs[i]-k,  k++) 
Fsc[pos4]  /«  sqrt2; 

} 

/♦ 

*  solve  least-squeures  problem;  need  workspace  of  size  m  +  nb*sz 

*  -  rhs  is  overwritten  by  dx 

*  -  in  first  Iteration,  estimate  condition  number  of  Fsc 
♦/ 

dgels„("N”,  &SZ,  &m,  ftintl,  Fsc,  ftsz,  rhs,  ftsz,  temp,  ftltemp, 

&inf o2) ; 


Fig.  4.  A  segment  of  C  code. 


for  (  int  i  =  0,  pos=0;  i  <  j;  pos  +=  j-i,  i++) 

{ 

double  seal  =  sigx(i); 

dx(pos)=  (  1.0/scal  +rho*scal)  *  sq2; 

} 

/*  loop  over  Fi,  coiiq)Ute  V’  *  Fi  ♦  V,  store  it  in  Fsc  */ 
for  (  vector  <  LaVectorDouble>: :  iterator  i  =  Fi->begin()+l; 
i  <  Fl->end() ;  i++,  n++) 

{ 

LaVectorDouble  to^  (Fsc .  addr  (}+pd.8z*n,  1 ,  pd^sz)  ; 
DualSccLle  (0 ,  *i,  vecx,  tmp)  ; 

/*  correct  diagonal  elements  ♦/ 
for  (  int  k  «  0,  pos*0;  k  <  j;  pos  +*  j-k,  k++) 
tnqp(po8)  *=  8q2; 

} 

LaLeast Square (dx,  Fsc,  ftn); 


Fig.  5.  A  segment  of  C++  code 

algorithm  for  semidefinite  programming.  The  program  also  uses  the  iterator  object  in  STL 
(Standard  Template  Library)  [10]  [8]  to  traverse  arrays  of  objects. 

The  first  major  difference  between  the  C  and  C++  implementations  is  the  way  the  initial 
data  is  read  in.  Unlike  C/Fortran  style  subroutines,  in  which  one  can  pass  a  pointer/address 
for  a  piece  of  storage  and  let  the  subroutine  split  the  storage  into  pieces  to  get  the  data, 
C++  objects’  constructors  have  no  such  scheme.  Initialization  is  done  by  reading  a  data 
file. 


An  Object-oriented  Approach  to  Semidefinite  Programming 


9 


The  second  major  diflFerence  is  that  because  C-H+’s  objects  are  higher  level  abstractions, 
the  implementation  in  C++  is  less  dependent  upon  pointer  arithmetic,  as  shown  by  the 
code  segments  in  C  and  C++  (Figs.  4  and  5)  for  doing  the  same  computation.  There  is 
overhead  associated  with  this  higher  level  of  abstraction,  but  we  will  show  that  the  effect 
on  performance  is  negligible. 

5  Conaparison  and  discussion  of  results 

Two  sets  of  data  are  obtained  by  randomly  generating  all  the  matrices  Fq,  Fi,  -  •  *,  Fm, 
and  the  vector  c.  Strictly  feasible  initial  points  so  and  Zo  are  also  generated.  The  timing 
results  are  shown  in  Table  1.  All  the  timings  are  done  on  a  HP  712/60  workstation.  Both 
the  C  implementation  from  [13]  and  the  C++  implementation  are  compiled  using  the  Gnu 
C/C++  compiler  version  2.7.2  with  the  same  compiler  options.  It  is  clear  from  Table  1 
that  the  performance  penalty  for  using  C++  is  only  a  few  percent  and  decreases  as  the 
problem  size  increases. 


Table  1.  Comparison  of  implementations. 


Example  1, 

n  =  40 

C++  implementation 

C  implementation 

time  (sec) 

time  (sec) 

20 

4.7 

4.4 

30 

6.9 

6.5 

Example  2, 

n  =  100 

50 

229 

222 

75 

390 

381 

We  have  shown  an  objected-oriented  design  and  implementation  of  a  semidefinite  pro¬ 
gramming  algorithm.  Even  though  object-oriented  technology  is  being  used  more  and  more 
widely  in  industry  now,  there  are  not  many  realistic  applications  to  numerical  computa¬ 
tion.  The  programming  environments  and  tools  seem  to  be  mature  enough  to  apply  this 
new  methodology,  and  the  performance  seems  to  be  comparable  to  a  non-object-oriented 
implementation. 

However,  there  are  other  considerations  that  have  to  be  taken  into  account  when  apply¬ 
ing  object-oriented  technology.  First,  it  takes  time  and  effort  to  learn  the  new  methodology. 
Second,  it  is  not  a  trivial  task  to  set  up  the  environment:  compiling  all  the  C++  packages, 
and  verifying  that  they  work  correctly,  especially  when  most  of  the  C++  packages  for  nu¬ 
merical  computation  are  still  in  the  testing  stage.  Third,  the  resulting  code  size,  i.e.,  the 
size  of  the  C++  executable,  is  about  2.5  times  that  of  the  C  executable.  With  continuing 
development  of  object-oriented  technology  and  of  compilers  for  object-oriented  languages, 
the  second  problem,  will  likely  be  alleviated.  The  hardware  considerations  of  the  third 
problem  are  becoming  less  of  a  hindrance  with  the  advances  in  the  computer  industry.  We 
do  believe  that  the  benefits  of  using  object-oriented  methodology  outweight  the  currently 
existing  disadvantages. 

The  design  and  analysis  in  this  paper  can  be  generalized  to  apply  to  object-oriented 
design  and  implementation  of  other  interior  point  algorithms  which  use  potential  reduction 
to  find  the  optimum. 
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